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ABSTRACT 


Described is the behavior of a rotorcraft equipped with Higher Harmonic Stall 
Control (HHSC) and a Reverse Velocity Rotor (RVR). Current rotorcraft are limited in 
forward flight speed by retreating blade stall and compressibility effects on the advancing 
blade. Stall occurs as the blade encounters increasingly severe reverse flow. HHSC 
enables conventional rotor systems to fly on the forward and aft sections of the rotor disk, 
greatly reducing reliance on the mixed flow regions defined by the advancing and 
retreating blades. Employment of the RVR allows lift generation while the rotor is 
experiencing reverse flow. A similar type of two per revolution (2/rev) input can be 
tailored to deliver maximum benefit to RVR equipped rotorcraft. 

Modification of the Joint Army Navy Rotorcraft Analysis and Design (JANRAD) 
computer program allows 2/rev cyclic input, use of the RVR, and analysis using high 
fidelity graphical output to examine angle of attack, coefficient of lift, and air load. 
Computational results show performance gains in conventional helicopters and high 
speed flight potential for RVR equipped aircraft. The RVR is applied to the Joint Heavy 
Vertical Lift (JVHL) aircraft conceptual design for preliminary analysis. This conceptual 
design can be used as an indicator of the performance of a high speed RVR equipped 
aircraft. 
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I. 


INTRODUCTION 


A. MOTIVATION FOR RESEARCH 

The next generation of ship-based military vertical heavy lift rotorcraft must be 
capable of conducting high speed long-range combat assaults. It must be able to operate 
in all weather conditions, day or night, from prepared or unprepared surfaces. 
Operational requirements are shown in the figures below. 

Long Range Heavy Lift Aircraft Mission Profile 



Figure 1. Mission Profile. 


Vehicle Gross Weights for Air Transport 
70 I-,-^^^-,- 



HMMWVLW155 M198 M7VR LAV HEMAT AAAV M1A1 
vehicle 


Figure 2. Payload Requirements. 
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These requirements are outside the envelope of any existing Vertical/Short Take- 

Off and Landing (VSTOL) aircraft to date. By the year 2020, the CH-46E will have been 

replaced by the MV-22A. But, the MV-22A does not have sufficient internal or external 

payload capacity to transport heavy equipment or large quantities of bulk logistics 

support. The CH-53E does not have sufficient internal load capacity to carry heavy 

equipment such as the Mk28/Mk928 7-ton truck needed to support USMC field artillery 

ashore. It does not have sufficient range with large external loads like the Eight Armored 

Vehicle (EAV) to support the Operational Maneuver from the Sea (OMETS) and Ship to 

Objective Maneuver (STOM) operations up to 200 NM inland envisioned in the future. 

Although a modified CH-53X has been proposed, no manufacturing plan has been 

formalized and current plans call for the phase out of existing CH-53Es by 2025. 

Clearly, a new aircraft must be developed to meet future requirements. 

B. HIGHER HARMONIC STALL CONTROL AND REVERSE VELOCITY 
ROTORS 

Two of the most compelling technologies that might be used to meet these 
operational requirements are Higher Harmonic Stall Control (HHSC) and Reverse 
Velocity Rotors (RVR). Both of these technologies have been studied in great detail but 
have never enjoyed widespread acceptance within the helicopter industry. An aircraft 
using HHSC and the RVR system would retain all the characteristics of a conventional 
helicopter at low airspeeds but would also be capable of high speed flight (in excess of 
200 knots). 

HHSC differs from Higher Harmonic Control (HHC) in the following manner. 
HHC is an active control concept whose principal objective has been to reduce helicopter 
airframe vibrations. More recently it has been applied to improve helicopter performance 
as well. Eor HHC principle frequencies of interest are (n-l)Q , nCl , and (n-l-l)Q in the 
rotating system and nO. in the fixed fuselage system. Here, n , is the number of blades 
and Q is the rotor rotational speed. [Wood et. al, 1985] In contrast to HHC, HHSC 
relates usually to 2fl or second harmonic frequency. [Arcidiacono, 1961] More recently 
it has been extended to 3Q frequency as well. The objective of HHSC when applied to 
the rotor in high-speed flight is to smooth out the distribution of lift and eliminate stall 
over the entire rotor disk. 
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A preliminary design of an RVR based vertical heavy lift aircraft was developed 
but analysis and verification of rotor performance was nearly impossible since traditional 
design methodology and performance calculations worked well with the rotorcraft only 
up to 170 knots. The planned cruise speed of the RVR based aircraft was 300 knots. 
After several attempts to constrain our design to work with traditional rotor design tools it 
was determined the best course of action would be to develop a new way to investigate 
the behavior of the rotor system in which both RVR and HHSC technology would be 
applied. 

C. JANRAD 2P_RVR 

Due to the unique nature of RVR and HHSC the original Joint Army Navy 
Rotorcraft Analysis and Design (JANRAD) program was modified to allow for HHSC, 
use of RVRs, and several other enhancements. The core pitch-thrust moment iteration 
technology that forms the core of JANRAD was unchanged. 



Figure 3. Graphical Depiction of Pitch-Thrust Moment Iteration. [From: 

Gerstenberger and Wood, 1963] 

A complete description of the original program, assumptions, and operating 
constraints can be found in theses by Nicholson [Nicholson, 1993] and Eccles. [Eccles, 
1995] Sections two and three discuss the employment and effects of using the HHSC and 
RVR systems. These chapters will provide a working knowledge of how each system 
works by outlining program modifications, assumptions, and actual JANRAD 2P_RVR 
output. UH-60 aircraft data will be used to further illustrate applicability to current 
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helicopter design. JANRAD 2P_RVR was used for the analysis of the Joint Vertical 
Heavy Lift Aircraft (JVHL) outlined in section four. Complete user instructions for 
JANRAD 2P_RVR are provided in section five. 

The study concludes with a brief list of recommendations for improvement and 
suggestions for further study. A complete listing of all the modules required to run the 
JANRAD 2P_RVR program is given in the appendices. 
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II. HIGHER HARMONIC STALL CONTROL 


A. INTRODUCTION 

Since the early 1960s researchers have pursued the idea of incorporating two per 
revolution (2/rev or 2P) control inputs to rotor systems. Studies ranged from increasing 
forward fight speed to the reduction of aircraft vibrations. Computer controlled Higher 
Harmonic Control (HHC), 3P, 4P and 5P, showed significant vibration reduction and was 
the subject of a detailed flight test program conducted by Hughes Helicopters, Inc. under 
a NASA/Army contract. [Wood et. al, 1985] For the study of this thesis the emphasis 
will be on performance airspeed increases and controlling the rotor stall pattern in which 
2P and 3P Higher Harmonic Stall Control (HHSC) will be applied. 

B. BLADE PITCH MOTION AND THE SWASHPLATE 

Blade pitch motion comes from two sources, namely commanded collective and 
cyclic pitch input from the helicopter control system and elastic deformations (twist) of 
the blade and control system. This section will focus on the cyclic commanded input. 
Elastic deformations are neglected at this level of study. 

The pitch time history of the rotor blade can be represented by a Fourier series as 
follows: 

6 iy/) = 6 q + ^ 1 , cos(^^) + sin(^) +... + cos(m//) + sin(n^^) +... 

where 0 is the blade’s feathering angle, y/ is the blade azimuth position, and y/ = Q.t. 
Conventional control inputs produced by the swashplate consist of collective pitch, 0 ^, 
and the first harmonics of the Fourier series: the lateral cyclic 0^^ and the longitudinal 
cyclic A conventional blade pitch (feathering) motion can be described by the 
Fourier series 

0 (y/) = 0Q + 6 * 1 ^ cos(^i^) + sin(^) = 6*0 + 0^^ cos(D.t) + 0^^ sin(Q 0 

This one per revolution control (1/rev or IP) formed the core of JANRAD’s 
original trim routine. The 0^^ term controls the lateral orientation of the rotor disk and 

6 *;^ term controls the longitudinal orientation. In the above equation, t = 0 has been 
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taken with the blade in the trailing or aft {y/ = Q) position. Remember that for a rotor 
with a centrally located flap hinge there is an exact 90 degree force/displacement phase 
lag. In the case of pure 6^^ (sine) longitudinal cyclic motion, the minimum aerodynamic 
flapping motion occurs at ^ = 90 degrees, where the minimum flapping displacement 
occurs 90 degrees later at ^ = 180 degrees; the maximum aerodynamic flapping moment 
occurs at y/ = 270 degrees, and the maximum flapping displacement occurs 90 degrees 
later, at = 360 degrees. 

The swashplate is the key to affecting the pitch control of the rotor blades. The 
swashplate has rotating and nonrotating disks concentric with the rotor shaft. Both disks 
slide up and down in response to collective inputs or they can be tilted to an arbitrary 
orientation in response to cyclic inputs made by the pilot. The pilot cyclically changes 
the pitch of the blades about the swashplate. HHSC inputs can be incorporated into 
existing helicopters by modification of the swashplate. HHSC inputs would be 
transmitted to the rotor using a stationary outer member with a track attached to a 
conventional swashplate assembly. This differs dramatically from HHC systems which 
rely on computer controlled actuators connected to the swashplate to alter the rotor pitch. 
Figure 4 compares a HHC system with a proposed HHSC mechanism. 



Figure 4. HHC Actuator System [From: Wood et. al, 1985] Compared to a HHSC 
Swashplate and Ring System. [From: Arcidiacono, 1961] 
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C. EXPERIMENTAL ANALYSIS USING JANRAD 2P_RVR 

Since research by Arcidiacono has suggested a significant increase in maximum 
forward speed of a helicopter could be obtained through the use of HHSC control of rotor 
pitch the first modifications to JANRAD were made to the Fourier series that defines 
cyclic pitch. First iterations of JANRAD 2P_RVR used the following equation for 
0 given by 

6 iy/) = 6 q+ 6'^ cos(^/^) + 6'^ sin(^/^) + 6^^ co%{2y/) + 6^^ %m{2y/) + sin(3^/^ + 

This equation incorporates lateral and longitudinal 2P signals and longitudinal 3P 
signals plus a phase, . The 2P and 3P 6 coefficients are user defined fixed values 

while the IP coefficients are determined by JANRAD based on first harmonic trimming 
rules. Since the objective of our study was to maxi miz e performance it was found that 
the best equation for 6 was found to be 

0 iy/) = 6»o + cos(^^) + sin(^^) + 6^^ cos(2ij/) + 6^^ sin(3^^ + ) 

The 2P sine term was eliminated because it was causing increased rotor stall on 
the left and right side of the rotor disk, in the exact location were we wished to reduce 
rotor stall. The 3P sine term produced striking results and was retained for experimental 
value, but in the majority of cases this term was assigned a value of zero. 

D. TEST USING PROUTY’S HELICOPTER 

As stated earlier, the primary reason for using HHSC (2P) is to improve the lifting 
efficiency of the rotor by modifying the stall pattern. Before attempting to analyze 2P 
behavior with JANRAD 2P_RVR, conventional first harmonic inputs were verified by 
comparing basic output with known results. Prouty’s sample helicopter was used as the 
baseline for comparison. To determine if JANRAD 2P_RVR was producing correct 
output, rotor disk angle of attack contour plots were compared to known solutions. These 
plots were chosen because the stall patterns are directly linked to the angle of attack 
distribution. All 2P and 3P terms were set to zero for the test. Figure 5 shows the 
calculated first harmonic cyclic pitch of the Prouty example helicopter. 


7 




Figure 5. Rotor Azimuth versus Harmonic Cyclic Pitch for the Prouty Example 

Helicopter. 

Figure 6 compares the angle of attack contour plot from Prouty’s text to JANRAD 
1P2P3P output. 


« - iw 



Figure 6. Prouty’s Example Helicopter Contour Plot (Feft) Compared to JANRAD 

2P_RVR Output (Right) Using Prouty’s Example Helicopter. 

The similarities between the two plots indicate JANRAD 2P_RVR output is 
stable and in close agreement with Prouty’s numerical results. Next, 2P and 3P signals 
were introduced. Figures 7 and 8 show the result for a 2P signal and 3P signal, 
respectively. 
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Figure 7. 2P Signal Effects on Cyclic Pitch (Left) and Angle of Attack (Right). 



Figure 8. 3P Signal Effects on Cyclic Pitch (Left) and Angle of Attack (Right). 


Three dimensional lift distribution plots were used to further verify the effects of 
these higher harmonic signals. The plots are compared to conventional lift distributions 
in Figures 9 and 10. 



Figure 9. Conventional (IP) Lift Distribution (Left) Compared to 2P Lift 

Distribution. 
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These plots clearly show how the 2P signal moves the aerodynamic forces to the 
forward and aft sections of the rotor disk. By doing this the helicopter can ‘fly’ on the 
forward and aft sections of the disk and therefore not rely as heavily on the left and right 
sections. The rotor is no longer trying to provide propulsive and lifting force from the 
unstable left and right rotor regions. This translates into the ability to gain forward 
velocity. As was mentioned earlier, the 3P signal was ineffective in showing any 
performance gains but was retained for experimental value. 

E. COMPARISON TO ARCIDIACONO 

One additional test was run to verify JANRAD 2P_RVR’s integrity. Thrust 
distribution and angle and angle of attack distribution were compared to ideal thrust and 
angle of attack distributions derived by Arcidiacono. [Arcidiacono, 1961] Figures 11 and 
13 show Arcidiacono’s original plots and his attempts to meet the optimum distribution. 
Note that Arcidiacono’s sign conventional is opposite to the sign conventional used 
throughout this study. 
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THRUST COEFFICIENT, ONE BLADE 



AZIMUTH ANGLE , , 0E6 


/I >.34 



Figure 11. Arcidiacono’s Results with No 2P. 


115.1492 Kts. advance ratiQ=0.29772. 2P input=0 degs. TPP angle=0.090769 degs. 3P input=0 degs. 3P offset=0 degs 


Rotor Tip Path Plane Stall Pattern, advance ratio=0.29772 



X 'lO'^ Actual and Ideal Thrust Distribution at Stall, advance ratio=0.29772 



Figure 12. Results Using JANRAD 2P_RVR with No 2P. 
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0 9 0 180 270 360 
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Figure 13. Arcidiacono’s Results with 2P and 3P. [From: Arcidiacono, 1961] 


115-1492 Kts. advance ratio=0-29776- 2P inpLit=4 degs. TPP angle=0.089514 degs. 3P input=0 degs. 3P offeet=0 degs 
Rotor Tip Path Plarte Stall Pattern, advance ratio=0.29776 



X 10"^ Actual and Ideal Thrust Distribution at Stall, advance ratio=0.29776 



Figure 14. Results Using JANRAD 2P_RVR with 2P. 


Figure 14 shows that JANRAD results were favorable in terms of thrust 
coefficient. Angle of attack distribution was expected to be less favorable due to the 
constraints imposed by JANRAD for tip losses. Overall, JANRAD 2P_RVR validated 
itself for use with 2P inputs. 
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F. FORWARD FLIGHT VELOCITY INCREASES 

Previous versions of JANRAD had been validated for use with the UH-60A; 
therefore it was selected as the helicopter to test for potential speed increases. [Eccles, 
1995] All UH-60A parameters were loaded into JANRAD 2P_RVR. Sikorsky’s 
proprietary airfoil, the SC1095R8, was used as the rotor blade. Testing was confined to 
conventional helicopter configurations only. 

In short, an increase in maximum flight speed of 14 knots was achieved using a 
four degree fixed 2P input. This is just short of a ten percent improvement. No 
significant difference in maximum forward speed was observed with 2P inputs greater 
than 4 degrees. This increase was observed under the 1.5 percent iterative deviation rule 
imposed by JANRAD 2P_RVR. That is, collective and cyclic pitch settings are adjusted 
until the thrust vector is within 1.5 percent of the previous iteration. This stringent 
requirement ensures realistic simulation within the trim module. 

The maximum airspeed without any 2P input was 159 knots. The rotor is under a 
great deal of stress to produce all the required forces for flight. Figures 15 through 18 
depict the rotor behavior under these conditions. Any increase in speed caused JANRAD 
2P_RVR to diverge and return to the user a ‘WILL NOT TRIM’ message. 
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Figure 16. UH-60A Lift Distribution at 159 Knots. 
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Figure 17. UH-60A Angle of Attack Distribution at 159 Knots. 
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Figure 18. UH-60A Cl Distribution at 159 Knots. 

The maximum airspeed with the four degree 2P input was 173 knots. The 
familiar 2P pattern is clearly evident in Figure 21. The other figures show results similar 
in appearance to the plots derived from Prouty’s example helicopter. Any increase in 
speed caused JANRAD 2P_RVR to diverge and return to the user a ‘WILL NOT TRIM’ 
message. 
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Figure 19. UH-60A 2P Cyclic Pitch at 173 Knots. 


Lift Distribution At 173.2245 Kts, advance ratio=0.38913, 2P input=4 degs, 3P input=0 degs, 3P offsBt=0 degs 
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Figure 20. UFI-60A 2P Lift Distribution at 173 Knots. 
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Figure 21. UH-60A 2P Angle of Attack Distribution at 173 Knots. 
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CL Distribution At 173.2245 Kts. advance ratio=0.36913. 2P input=4 degs. TPP angle=19 1146 degs. 3P input=0 degs. 3P offset=0 degs 


Figure 22. UH-60A 2P Cl Distribution at 173 Knots. 


Although not shown, testing was also carried out with compound helicopters. Use 
of 2P systems with compound helicopters equipped with auxiliary thrust devices 
performed very similarly to conventional 2P equipped helicopters while the rotor is 
providing forward thrusting power. As the auxiliary thrust system begins to provide the 
majority of propulsive force at airspeeds above 180 knots, if no power was supplied to 
the rotor, the disk would assume an aft tilting tip path plane and rotor drag would 
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increase significantly. This situation is avoided completely in actual compound 
helicopters in high-speed flight by providing a relatively small amount of power to 
stabilize the rotor, while keeping the tip path plane at a near zero angle of attack. 

The application of second harmonic control is not a new idea for improving the 
speed performance of helicopters. This thesis has verified with JANRAD 2P_RVR that a 
2P rotor system does trim at airspeeds up to 10% higher than conventional IP rotor 
systems. The next section explores the use of Reverse Velocity Rotors coupled with a 
HHSC control signal. 
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III. EMPLOYMENT OE REVERSE VELOCITY ROTORS 


Current helicopters are limited in their speed by three effects, namely stalling of 
the tip of the retreating blade, loss of lift due to reverse flow over the inboard portion of 
the retreating blade, and compressibility on the tip of the advancing blade. [Ewans, 
1973] The severity of these effects is a function of the rotational speed of the rotor 
blades and the forward velocity of the aircraft. Note that if we treat retreating blade stall 
by slowing the rotor we also reduce compressibility effect on the advancing blade. 
Several studies have been conducted to deter mi ne the best way to combat retreating blade 
stall. Everything from exotic blade shapes, segmented rotors [Zientek, 2001] and forced 
air/circulation control have been attempted. 

A. REVERSE VELOCITY ROTORS 

A relatively new and promising concept recently introduced is called the Reverse 
Velocity Rotor (RVR). RVR systems are a recent development, invented by Harold 
Lemont to address the aerodynamic issues contributing to retreating blade stall and 
advancing blade compressibility phenomena, the primary factors that currently limit 
helicopter maximum forward flight speeds. Instead of eliminating the reverse flow that 
occurs on the retreating side of the blade the RVR is actually designed to operate in the 
reverse flow region. With conventional rotor blades, this would completely eliminate the 
potential to generate lift. In fact, a ‘reversible’ airfoil section with a rounded trailing edge 
is employed to give reasonable lift and drag characteristics in the reverse flow. Reverse 
velocity rotor systems are built around double-ended airfoils as shown in Eigure 23. 


Typical RVR airfoil 



Eigure 23. RVR Cross Section. [Erom: Ashby, 2002] 

Reverse velocity rotor systems represent a revolutionary approach to high-speed 
Vertical Take Off and Landing (VTOL) configurations. Application of this concept 
makes possible rotorcraft designs that are capable of attaining speeds significantly greater 
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than conventional or even compound helicopters. Initial analysis shows cruise speeds in 
excess of 300 knots are possible with RVR systems. Sikorsky Aircraft Corporation has 
completed a study that shows future Runway Independent Aircraft (RIA) equipped with 
RVR systems could attain cruise speeds over 300 knots, double the cruise speed of 
conventional helicopters. [Ashby, 2002] This airfoil design minimizes the impact of 
retreating blade apparent velocity reduction caused by summing the blades rotational 
velocity and the aircraft’s forward flight velocity as illustrated in Figure 24 below. 


Vapparent = 275 + 507 = 782 ft/S 



Figure 24. Rotor Disk Apparent Velocities. [From: Ashby, 2002] 

As might be expected the RVR system does require somewhat more power to 
hover and fly at low airspeeds. For this reason the best application for RVR systems is 
for rotorcraft that will spend the majority of time in cruise configuration, i.e. civilian or 
military transport. Aircraft that will spend the majority of their flight time in low speed 
flight or hovering, i.e. rescue aircraft, should continue to employ conventional rotor 
blades. Figure 25 shows predicted power requirements. 
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Figure 25. Power Requirements versus Velocity. [From: Zientek, 2001] 


B. WIND TUNNEL TEST 

Fairchild Republic Division conducted wind tunnel tests at NASA Ames in 1973 
of a 1/7 scale 12% thickness double ended airfoil at equivalent airspeeds of over 310 kts. 
[Ewans, 1973] The test rig is shown in Figure 26. 



Figure 26. Installation of RVR Test Rig in the NASA Ames 12 Foot Wind Tunnel. 

[From: Ashby, 2002] 
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Testing was considered a success although some large unbalanced forces were 
noted at the edge of the test envelope. Data collection focused on measuring rotor 
performance, particularly maximum lift and effective lift/drag ratio over the full 
envelope. [Ewans, 1975] The most important data for this study are given in Figures 27 
and 28. From these plots data tables were generated which correlated Mach number, 
angle of attack, and Cl or Cd. These tables were inserted into a MATFAB function 
which used a two-dimensional interpolation to provide the JANRAD 2P_RVR thrust 
moment and drag moment routines with the appropriate Cl and Cd values. 


sEcnoH urt cocrnciENT through M<r anglc of attack 
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Figure 27. RVR Cl versus Angle of Attack. [From: Ewans, 1975] 
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Figure 28. RVR Cd versus Angle of Attack. [From: Ewans, 1975] 

C. EXPERIMENTAL ANALYSIS USING JANRAD 2P_RVR 

JANRAD, in its original form, was not set up to calculate trim parameters for a 
rotorcraft traveling at 300 knots. Several new assumptions had to be made to make 
JANRAD 2P_RVR compatible with high speed RVR systems. 

I. Advance Ratio 

Before outlining the assumptions, the advance ratio parameter must be defined. 
The advance ratio is the ratio of the rotor’s horizontal speed to its tip speed. The tip 
speed ratio, /u , is defined as: 

A = Kd cos a / QR 

where 

= forward airspeed 
a = tip path plane angle 
Q = rotational velocity in radians per second 
R = rotor radius 

With a conventional rotor system as ju exceeds 0.4, collective and cyclic pitch 
values were unable to satisfy rotor lift, rotor drag, or resultant thrust vector requirements. 
The system falls out of ‘trim’ with respect to pitching, rolling and yawing moments. Even 
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after applying second or HHSC the program was unable to trim above advance ratios of 
.5. This is not surprising when we recognize that at // = 0.5, one half of the retreating 
blade {y/ = 270 degrees) is in reverse flow. In contrast, RVR systems operate in an 
environment with advance ratios ranging between 1.0 and 2.0. This fivefold increase in 
advance ratio in RVR aircraft over aircraft equipped with conventional rotors illustrates 
the significant differences between the two systems. 

2. RVR Assumptions and Requirements 

Several assumptions were made to enable computational trimming of an RVR 
system using the traditional JANRAD trimming program. For high-speed rotorcraft 
flight (// > 0.5) JANRAD 2P_RVR only allows RVR trimming on compound helicopters 
where a wing provides the majority of vertical lift at cruise speeds. Proper employment 
of a wing can reduce rotor loads up to 80 percent. The reduction in rotor loading 
eliminates blade stalling concerns due to high loading. It is also assumed that the rotor of 
a RVR helicopter will operate in or near to an auto-rotative condition (small amount of 
rotor power at low tip path plane angle) with auxiliary propulsion of the vehicle. [Ewans, 
1973] While in this flight condition the stability axis of the aircraft is moved away from 
the center of the rotor disk to a location closer to that of a typical fixed wing aircraft. 
This relieves the lateral and longitudinal cyclic pitch constraints since the rotor disk will 
seek an equilibrium condition based primarily on collective pitch. The final condition 
that did not allow the RVR system to trim at high speed was rotor drag. With the RVR 
system in full operation (// = 2 and negative pitch on the retreating side of the rotor disk) 
conventional helicopter rotor drag analysis does not work. Fortunately, having the rotor 
in a near auto-rotative state allows a simple but critical drag assumption; over the entire 
rotor, the integrated effect of the tilt of the lift vector at each blade element is sufficient to 
overcome the integrated effect of the drag at all of the blade elements. [Prouty, 1986] 
Only a small amount of power is required to overcome the drag and provide rotor 
stability. These assumptions are highly interdependent. All of the conditions must be 
satisfied concurrently whenever RVR trim calculations are applied. 

Since having the correct advance ratio, //, is critical to successful RVR operation 
JANRAD 2P_RVR has a built in calculator which computes the optimal // for the design 
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cruise airspeed such that the entire retreating blade is in reverse flow. This optimal ju 
value ensures the reverse flow area around the retreating blade is sufficient to support 
RVR aerodynamic operation and RVR based cyclic pitch, discussed below. 

3. RVR Based Cyclic Pitch 

One final element had to be added to the RVR system: RVR based cyclic pitch. 
The use of HHSC pitch is required for advance ratios greater than 1.0. [Zientek, 2001] 
Using 1/rev cyclic pitch alone, the rotor cannot achieve the required large negative pitch 
in the 270 degrees blade sector, yet maintain positive pitch everywhere else. See Figure 
29. 


loHow Angle (deg) 



Figure 29. Inflow Angle versus. Blade Station. [From: Zientek, 2001] 

The 2/rev RVR based cyclic pitch can obtain the proper blade angle of attack on 
the retreating side without the aft tilt normally associated with windmilling operation. 
This differs from conventional HHSC. The equation for 6 takes the form 

= 6»o + 6»o sin(^^) - 0 ,^ cos(^^) - 0 ,^ sin(^^) + 0 ^^ cos(2^^) 

The negative signs preceding the 6^^ and 6^^ terms are required to move the 

maximum cyclical inputs to the retreating sections of the rotor disk. Figure 30 compares 
conventional HHSC based cyclic pitch and RVR based cyclic pitch. 
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Figure 30. Conventional HHSC Cyclic Pitch (Left) Compared to RVR Based Cyclic 

Pitch. 


D. EXPERIMENTAL RESULTS 

All of these assumptions are enforced within JANRAD 2P_RVR by implementing 
user enabled logic gates. The logic gates allow JANRAD 2P_RVR to operate in either 
conventional or RVR based cyclic pitch. JANRAD 2P_RVR only contains one non-user 
driven restriction for RVR based cyclic pitch employment, airspeed must be greater than 
200 knots. Figures 31 through 35 show JANRAD 2P_RVR output. 

Figure 31 clearly shows the effects of the RVR Based Cyclic Pitch. The negative 
blade angles at rotor azimuth positions between 230 degrees and 320 degrees are crucial 
to the development of positive lift. Also note that collective pitch has been reduced to 
zero through JANRAD 2P_RVR’s trim constraints, complementing the assumption that 
the rotor system has entered a near autorotative state. 
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Figure 31. Nature of RVR Based Cyclic Pitch. 


Figure 32 shows the Mach number distribution around the rotor disk. Note that 
the retreating side of the blade, approximate blade azimuth positions of 200 degrees to 
340 degrees, is fully entrenched in reverse flow. It is the combination of this reverse 
flow, the RVR’s local angle of attack, and Cl that allow for positive lift distribution. 
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Figure 33 shows Cl distribution viewed as a contour plot. One can immediately 
see that the disk is trying to achieve equilibrium. FIHSC contour lines are replaced with 
the RVR based cyclic pitch contour lines. The RVR based cyclic pitch builds upon the 
principles of FIFISC but shifts the emphasis on stall control from the forward and aft 
portions of the disk to the left and right sides of the disk. The shift of stall control to the 
lateral sections of the disk is critical to take full advantage of the traits of the RVR. The 
contour lines in Figure 33 clearly illustrate that the RVR based cyclic pitch is functioning 
as predicted. Also apparent in Figure 33 are the high Cl values in the forward right 
sector of the disk. This phenomenon is discussed in the following paragraphs. The 
distinct left and right side loading is further illustrated in Figures 34 and 35. 
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CL Distribution At 300.3893 Kts, advance ratio=2.0003. 2P input=9 degs. TPP angle=3 degs. 3P input=0 degs. 3P offset=0 degs 



90 


Figure 33. Cl Distribution Contour Plot. 


Figure 34 shows the radial distribution of blade lift at the most important rotor 
azimuth locations. At the y/ = 270 degree location positive lift is observed even though 
the blade is total reverse flow. As stated earlier, the positive air load is a result of the 
RVR’s unique aerodynamic behavior given the proper operating conditions. 
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Figure 34. Blade Position versus Air Load at Critical Rotor Azimuth Locations (0, 

90,180, and 270 Degrees). 


Figure 35 is a three dimensional representation of the lift distribution around the 
disk. Note the predominant lift ‘hump’ centered at the disk’s 270 degree azimuth 
location, further proof the RVR based cyclic pitch is performing as predicted. Also note 
the ‘spike’ centered at approximately ^ = 160 degrees. This ‘spike’ is the result of high 
blade angle, and therefore high blade coefficient of lift. The high blade angle is due to 
the cumulative effects of the cyclic inputs within that azimuth quadrant. Although 
undesirable, this phenomenon does not effect the disk’s ability to trim within the 
computational boundaries established in JANRAD 2P_RVR. 
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Figure 35. Lift Distribution Surface Plot. 


This area of abnormal lift distribution will certainly cause vibrations and high 
blade stress as each rotor blade is forced to pass through the effected azimuths. In fact, 
Ewans encountered similar unbalanced lift forces in his tests [Ewans, 1975] that made it 
impossible to complete his entire test program. In an effort to smooth the lift distribution 
a 3P signal (4 degrees) was added. The resulting lift distribution is shown in Eigure 36. 
In this figure three distinct humps are visible and the overall lift distribution is generally 
equalized between the three humps. Application of additional higher harmonic inputs 
would continue to suppress unwanted lift spikes. Detailed study of rotor system behavior 
using harmonics higher than 3P is beyond the scope of this study and would require 
further analysis in future research. 
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Lift Distribution At 300.3893 Kts, advance ratio=2-0003.2P input=3 degs, 3P input=4 degs, 3P offset=0 degs 



Figure 36. Lift Distribution Surface Plot with 3P Harmonic Cyclic Pitch. 

JANRAD 2P_RVR has shown that a RVR system with tailored 2P cyclic pitch 
does trim under the assumptions used in the study. The ability to trim allows for 
sustained high speed aerodynamic flight without incurring high drag penalties or autogiro 
type behaviors. The ability to maintain aerodynamic flight with a rotor system above 200 
knots allows designers to expand the performance and usability of current helicopters. 
The next section explores one possible way to employ a complete RVR system to fulfill 
an emerging vertical heavy lift requirement. 


32 
































IV. APPLICATION OF THE RVR AND ACCOMPANYING 
CONTROL SYSTEM TO THE JOINT VERTICAL HEAVY LIFT 

(JVHL) AIRCRAFT 

A. FIRST ITERATION DESIGNS TO SATISFY REQUIREMENTS 

In an effort to satisfy operational requirements defined in the first chapter, two 
‘new-start’ aircraft, the Joint Heavy Lift (JHL) and Reverse Velocity Rotor (RVR), were 
the subjects of simultaneous preliminary design analysis. For more information on these 
advanced designs the reader is referred to the references entitled “Joint Heavy Lift 
Expeditionary Warfare Aircraft Solution” [Wood and Aaron, 2003] and “Reverse 
Velocity Rotor Approach for Design of a STOVL Heavy Lift Transport for Expeditionary 
Warfare.” [Wood and Van Riper, 2003] The JHL is intended to be a long range, heavy 
lift aircraft to support Marine Corps, joint, and coalition force operations ashore up to 200 
NM inland in a forcible entry environment. The JHL is a Joint Strike Tighter (JSE) lift- 
fan variant modified for the increased load-carrying requirement. Large lift fans, each 
producing up to 60,000 pounds of vertical thrust, are embedded in the wing sections of 
the aircraft to provide Vertical Take-off and Landing (VTOL) capability while thrust 
vectored turbofan engines provide forward thrust while in forward flight. The RVR 
aircraft is a compound helicopter equipped with a conventional anti-torque tail rotor and 
large turboprop engines to provide the auxiliary thrust required for high-speed forward 
flight. The rotor system is an eight bladed (110 foot diameter), fully articulated, foldable 
system with each blade incorporating the RVR airfoil cross section. The rotor system is 
driven by a lightweight variable speed transmission that allows RPM to be set at either 
100% or 50% of normal operating RPM based on flight mode. The anti-torque system is 
a traditional six bladed pusher tail rotor mounted to the tail vertical pylon. Auxiliary 
propulsion is provided by two propellers (fifteen foot diameter), each driven by one of 
the two wing-mounted turboprop engines. Eigures 37 and 38 show each design. 
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Figure 37. RVR Based Aircraft. 



Figure 38. JFIL Aircraft. 


Both aircraft are intended to operate from ships of the amphibious task force, 
specifically the Amphibious Assault and Logistics Support variants of the family of 
ExWar platforms currently under design by the Total Ship Systems Engineering (TSSE) 
curriculum at the Naval Postgraduate School. The aircraft are also designed to be 
compatible with legacy force ships, i.e. LHD, LHA (R), and MPE (E) class ships to the 
maximum extent possible. 

As might be expected these evolutionary designs became the subject of much 
debate and criticism. Concerns about the RVR ranged from the expected size, weight, 
and complexity of a variable speed main transmission able to handle an aircraft with a 
take-off gross weight of over 100,000 pounds to the employment of a tail rotor almost as 
large as a UH-60 main rotor blade. The JFIL suffered from lift fan and wing integration 
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challenges and a downwash velocity, the vertical wind produced as the aircraft hovers, 
which could practically lift concrete blocks off the ground and send them flying several 
yards away. Clearly, a more amiable solution had to exist. 

B. THE JOINT VERTICAL HEAVY LIFT AIRCRAFT 

Ongoing research has produced a new conceptual design that combines the best 
traits of the RVR and JHL aircraft while greatly reducing the negative aspects of each 
design. Although the JHL did not employ any RVR or HHSC technologies that form the 
core of this study, it provides the framework for the FI 19 (JSF) Shaft Driven Lift Fan and 
Propulsion module integration. This technology provides the key to overcoming the most 
challenging design point of the aircraft; sustained hover and low speed flight in harsh 
environments. The proposed hybrid design, labeled the Joint Vertical Heavy Lift (JVHL) 
uses a RVR based rotor system and drive train with three F-35 V/STOL type Shaft 
Driven Lift Fan and Propulsion Modules. A direct jet anti-torque thruster is used in lieu 
of a conventional anti-torque tail rotor. JANRAD 2P_RVR was used exclusively for the 
JVHL’s rotor performance calculations. A detailed description of JANRAD 2P_RVR 
and user instructions are provided in Chapter V. The proposed configuration of the 
JVHL is shown in Figures 39-42. Dimensions shown in Figure 39 are in feet. 



Figure 39. JVHL Three View. 
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■j C-130J sized fuselage 


Figure 40. Top View of JVHL. 
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Figure 41. 


Perspective View of JVHL. 


Figure 42 shows the JVHL in hover configuration with 60,000 lbs of lift from the 
main rotor, 18,000 pounds from each fuselage mounted FI 19 module and 18,000 pounds 
from the tail mounted shaft driven lift fan, for a total of 114,000 lbs of vertical lift. 
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Figure 42. JVHL in Hover Configuration. 

Immediate benefits of using the FI 19 propulsion modules are clearly evident. 
Most important, the propulsion module offsets the hover and low speed performance of 
the RVR system by providing a great deal of vertical thrust. Current FI 19 engine 
configurations would provide ample power for all flight modes of the JVHL. In fact, the 
three FI 19 modules would provide such a surplus of power that each module could be 
down-rated by up to 30%. This modification would undoubtedly reduce fuel 
requirements and, more importantly, reduce engine fatigue. 

1. Fuselage Modules 

The fuselage mounted modules are installed so that the drive shafts leading to the 
main transmission have a straight path and direct jet exhausts do not have to be ducted 
around the fuselage. When in forward flight the two fuselage mounted modules provide 
ample thrust so the JVHL can achieve the design speed of 300 knots while at an altitude 
of 10,000 feet MSL. The large compound wing no longer has to accommodate turboprop 
engines or complicated shafting. The wing can now have a full complement of control 
surfaces and even be modified to carry extra fuel. Of course the wing will still have a 
negative effect on hover performance but the added vertical thrust provided by the direct 
jet thrusters and the aft lift fan help reduce the wing’s adverse effects. 

2. Tail Module 

The tail module eliminates all tail rotor components and replaces them with a self- 
contained highly reliable vertical lift augmentation and anti-torque solution while 
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significantly reducing drag in high-speed flight. The vectored thrust nozzle is collective 
and pedal coupled to ensure maximum anti-torque thrust and directional control is 
delivered proportional to inputs. While hovering, the module would also provide 
approximately 20,000 lbs vertical lift via the built in Lift Fan. This lift would help offset 
main rotor lift and power requirements and extend the CG range of the aircraft. 

3. Preliminary Design Data 

Another benefit of the main rotor, direct jet thruster, and left fan configuration is 
that the downwash will be less than half that of the JHL and the direct jet thruster output 
will be almost immediately cooled by the entrained flow circulating through the main 
rotor system. One of the most important advantages of using the FI 19 module is that 
now the JVHL will have a measurable power and performance margin that will enable 
the aircraft to hover with a full load on a High/Hot day (4000 ft MSL and 95 degrees F). 
Critical design data are shown in Tables 1 and 2. 




High Hot. 4000^t. 95degree F, rho 
is 0.00192 

Standard Dai 



Nofinal drive 

Normal drive 

Parameters 

Notation 

JVHL (hover) 

JVHL (fwd 
flight) 

JVHL (hover) 

JVHL (fwd 
flight) 

Weiqhts (lbs) 

Empty 

W, 

62865 

62865 

62875 

62875 

Payload 


37500 

37500 

37500 

37500 

Fuel Capacity 


20000 

20000 

20000 

20000 

Empty/Gross ratio 

W./GW 

0.52 

0.52 

0.52 

0.52 

Gross weight 

GW 

120365 

120365 

120375 

120375 

Rotor/compound winq proportion 

Lift from direct jet/aft 
fan 

LF 

50000 


50000 


% of lift for rotor 

SLr.I.r 

100% 

20% 

100% 

20% 

of lift for winq 


0% 

80% 

0% 

80% 

Rotor parameters 

Thrust (lbs) 

T = %L..i„*GW 

70365 

'24073 

70375 

24075 

Radius (ft) 

R 

40 

40 

40 

40 

Chord (ft) 

c 

3 

3 

3 

3 

No. of blades 

b 

8 

8 

8 

8 

FWD Speed (kts), 

Okts means hover 

V.-j 

0 

300 

0 

300 

Omeqa (RPM) 

Omeqa 

150 

60 

150 

60 ' 

Figure of merit 

FM 

0.8 

0.8 

0.8 

0.8 

Disc area (ft*2) 

A = pi*R''2 

5027.20 

5027.20 

5027.20 

5027.20 

Solidity 

Sigma = 
(b*c)/(pi*R) 

0.191 

0.191 

0.191 

0.191 


Table 1. Basic JVFIL Design Data. 
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Hi 9 k Hot, 4000ft, d5de9rec F, rko is 
0.00102 

Steederd De? 


Nora^l drive 

Nornel drive 

Parameters 

Notation 

JVHL (hover) 

JVHL (fvd 
flight) 

JVHL (hover) 

JVHL (fwd 
flight) 

Tip speed at 0 deg 
rotor azimuth 

Vllr.i = 

Omeqa*R 

628.40 

251.36 

628.40 

251.36 

Tip speed at 90 deg 
rotor azimuth 

Viip.ii = 

Omeqa'R • Vf-i 

628.40 

758.36 

628.40 

758.36 

Tip speed at 180 deg 
rotor azimuth 

V = 

Omeqa'R 

628.40 

251.36 

628.40 

251.36 

Tip speed at 270 deg 
rotor azimuth 

Vi:p.z7i = 

Omeqa'R - V«.j 

628.40 

-255.64 

628.40 

-255.64 

Mai Mach lip 

Mai Mi:, 

0.57 

0.68 

0.57 

0.68 

Disc loading fpsfl 

DL = Tf A 

14.00 

4.79 

14.00 

4.79 

Induced velociti f^ps1 

V: = 

sgrtrOL/2‘rho1 

60.37 

Not applicable 

54.25 

Not applicable 

Advance ratio 

mu = 

V*.j lOmeqa'R 

0.00 

w 

2.01 

0.00 

w 

2.01 

P:.j...jrHP] 

P:.j...j(HP) = 
T’V: 1 550 

7724.02 

Not applicable 

6941.94 

Not applicable 

1st appros: 

P..i..ifHPl 

..(HP) = 

P:.j...j f FM 

9655.03 

Not applicable 

8677.43 

Not applicable 

Thrust coed. 1 sigma 

Ct = 

T1 A'tho*(omega 
'R1''2'sigma 

0.088 

Not applicable 

0.078 

Not applicable 

Torgue coed. 1 sigma 

C* 1 sigma 
(from ProutI) 
teit page 90-92) 

0.0270 

Not applicable 

0.020 

Not applicable 

Torgue 

Q = 

786083 

Not applicable 

721183 

Not applicable 

2nd approi: 

P-i„ifHPl 

calculated from 

P = 

Q*omegaf550 
vhere Q is from 
Cqfsigma curve 
in teit page 90- 
92 

22450 

Not applicable 

20597 

Not applicable 

plus 33X 

Final P reguired 

29859 

Not applicable 

27394 

Not applicable 


Table 2. JVHL Rotor and Power Calculations. 

C. RVR ANALYSIS 

The JVHL uses two emerging technologies, the RVR system and the FI 19 
propulsion module. All other major components of the aircraft are based on fairly mature 
technologies and processes. Compound helicopter design is well understood, the CH- 
53E rated trans mi ssion and drive line technology needed to satisfy main rotor lift 
requirements has been refined by the helicopter industry, and use of composites in 
fuselage design practices is becoming more commonplace with each passing year. The 
only remaining design obstacle is the RVR rotor system. 

A detailed aeroelastic and dynamic analysis of the RVR system is beyond the 
scope of this study but a preliminary aerodynamic analysis was conducted to ensure the 
RVR system would be viable at high airspeeds. JANRAD version lP2P3Prvr was used 
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to check rotor system performance. As discussed earlier, correctly manipulating the 
advance ratio is critical since the RVR system relies on the entire retreating side of the 
rotor system being fully entrenched in reverse airflow. For the JVHL an advance ratio of 
2 is used at the aircraft’s cruise speed of 300 knots. To achieve this advance ratio the 
rotational speed of the rotor system is reduced from 150 rpm to 60 rpm. This high 
advance ratio causes sufficient reverse airflow on the retreating side of the rotor disk, 
trailing edge to leading edge, to produce lift. The RVR is coupled with a RVR unique 
two per revolution (2/rev) cyclic pitch signal. See Figures 43-45. 


Rotor Azimutii ^igle v. RV^ Harmonic Cyclic Pitch for 300 3893 Kts. advance ratIo=2 0003. 2P input=9 degs 



Rotor Azimuth Angle (degs) 


Figure 43. RVR Based Cyclic Pitch for the JVHL. 
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180 

Mach Number Distribution At 300.3893 Kts. advance ratio=2.0002 



Figure 44. Mach Number Distribution at Cruise Speed (300 kts),Retreating Side of 

Blade is in Reverse Flow. 


CL Distribution At 300 3893 Kts. advance ratio=2.0003. 2P input=9 degs. TPP angle=3 degs. 3P input=0 degs. 3P offset=0 degs 


270 


0 



90 

1 


Figure 45. Cl Distribution at Cruise Speed (300 knots) for the JVHL. 


Since the large compound wing has unloaded the rotor (at cruise speed the wing is 

carrying 80% of aircraft weight and the rotor is carrying the remaining 20%), the rotor 
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disk can operate in an autorotative state with small cyclic input and no wind milling. 
Without this wind milling behavior there is almost no drag penalty. These figures 
illustrate the behavior of this rotor system is directly in agreement with predicted RVR 
behavior. The aircraft is able to operate at high cruise speed without the typical negative 
effects of the rotor disk. 


42 



V. JANRAD USER INSTRUCTIONS 


A. GENERAL 

Program installation, guidance for input requirements, and rotor performance 
output are discussed in the subsequent sections. These instructions assume that the reader 
has a working knowledge of MATLAB release 12 and higher. 

Joint Army Navy Rotorcraft Analysis and Design (JANRAD) version lP2P3Prvr 
was developed to allow introduction of 2/rev cyclic inputs, 3/rev cyclic inputs, compound 
wing sizing, unconventional vertical lift augmentation, and employment of the Reverse 
Velocity Rotor (RVR) system and subsequent RVR specific 2P control inputs. This 
version also includes a comprehensive plotting capability that gives the user high fidelity 
2-D, 3-D, and contour plots to aid in rotor performance analysis. 

JANRAD version lP2P3Prvr was built using the original source code found in 
JANRAD version 1 [Nicholson, 1993]. This approach ensured the original trim routine 
that forms the heart of program was unaltered by subsequent undocumented updates. 

The program was developed and tested on a Pentium 4 based laptop running 
MATLAB release 13 under Windows XP Professional. The laptop was equipped with 
512 megabytes of RAM and a twenty gigabit fixed drive. 

B. INSTALLATION 

To install JANRAD 2P_RVR, create a subdirectory named JANRAD and install 
the following JANRAD 2P_RVR modules: janrad.m, trim.m, dmcalc.m, hh02clcd.m, 
ool2clcd.m, perf.m, plotter.m, polarl.m, rvrclcd.m, scl095r8clcd.m, thrcalc.m, tmcalc.m, 
and vrl2clcd.m. Listings of the modules can be found in subsequent appendices. 

Once this installation is complete, open MATLAB and type JANRAD (with the 
proper drive designation) in the Current Directory dialogue box. All of the modules must 
be installed in the select directory on the computer’s fixed drive for the program to 
function. The program will not work if the user attempts to run the modules from a 
remote drive, compact disk, or network. The program is now ready to run. 
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Figure 46. Title Page. 


At the File Selection page, enter a 1 if you want to edit a previously stored data 
file, otherwise enter a 2. If you chose to edit a file, you will be prompted on the next 
page to enter the file name without the extension. If JANRAD 2P_RVR is unable to find 
the specified file you will be prompted to try again. Ensure your spelling is correct and 
you are in the proper directory. Once the file is loaded the edit menu, as seen in Figure 
47, will be displayed. Check to make sure the menu’s title corresponds to the type of 
helicopter you are analyzing, either conventional or compound. 
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COHPOUMD HELICOPTER 

EDIT MENU 

1. 

pressure altitude 

2. 

temperature 

3. 

airspeed 

4. 

gross weight 

5. 

numlier of blades 

6. 

blade radius 

7. 

blade chord 

8. 

hinge offset 

9. 

blade grip length 

10. 

blade tTnist 

11. 

blade wieght 

12. 

number of blade elements 

13. 

rotational velocity 

14. 

if azimuth sectors 

15. 

lift curve slope 

IS. 

airfoil 

17. 

collective pitch 

10. 

flatplate area 

19. 

vert projected area 

20. 

horizontal tail area 

21. 

horizontal tail span 

22. 

horizontal tail CL 

23. 

horizontal tail CDo 

24. 

vertical tail area 

25. 

vertical tail span 

26. 

vertical tail CL 

27. 

vertical tail CDo 

28. 

2P input 

29. 

3P input 

30. 

auxiliary thrust 

31. 

Tiling area 

32. 

wing span 

33. 

Tiring CL 

34. 

wing CDo 

35. 

Tiring efficiency factor 



0. 

NO CHANGES 



Select the patameter to change 

= 1 



Figure 47. Compound Helicopter Edit Menu Page. 


Select the number of the parameter you wish to change. The current value for that 
parameter will be displayed followed by a prompt to enter the new value. If you chose 
not to change the value, press <ENTER> and the previous value will be maintained. 
Once a new value is entered or <ENTER> is pressed the edit menu will be displayed 
again. Choose the number of the next parameter to change or 0 to exit the programs 
direct edit mode. If you chose to create a new data file, prompts for individual 
parameters will be displayed. Enter the value of each parameter when prompted. 

Enter the values for the specified parameters regardless of editing or creating a 
data file in the following manner: 

1. PA - Enter the pressure altitude in feet 

2. temperature - Enter the temperature in degrees Eahrenheit 

3. airspeed - Enter forward airspeed in knots. JANRAD automatically 
converts airspeed to feet/second. 

4. GW - Enter aircraft gross weight in pounds. 

5. number of blades - Enter the number of blades used in the rotor system. 
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6. radius - Enter the rotor blade radius in feet. Measure the radius from the 
center of the rotor hub to the tip of the bade. 

7. chord - Enter the average rotor blade chord in feet. 

8. hinge offset - Enter the effective hinge offset in feet. 

9. blade grip - Enter the distance from the center of the rotor hub to the 

beginning of the aerodynamic portion of the rotor blade in feet. 

10. blade twist - Enter rotor blade twist in degrees. You can enter blade twist 
as either a positive or negative value. JANRAD uses the absolute value of 
twist. 

11. blade weight - Enter the weight of a single rotor blade in pounds. 

12. rotor blade elements - Enter the number of rotor blade radial elements. 
Multiples of 20 are recommended, i.e. 20, 40, 80. As the number of blade 
elements increases program execution is slowed. JANRAD adds one 
additional element to account for tip loss. 

13. rotational velocity - Enter the rotor rotational velocity in radians/second. 

14. azimuth sectors - Enter the number of rotor disk azimuth sectors. This 
number must match the number of rotor blade elements if JANRAD 
created plots are to be used. The number may be different if no plots are 
required. 

15. lift curve slope - Enter the average slope of the linear portion of the rotor 
blade’s airfoil’s lift curve, where the lift curve is Cl versus alpha. 

16. rotor blade airfoil - Eour airfoils are given for use with conventional 
helicopters and five airfoils are given for use with compound helicopters. 
Enter the number that corresponds to your desired airfoil. 

17. collective pitch - Enter the estimated collective pitch at 0.7 r/R in degrees. 
JANRAD automatically converts the vale to radians. 

18. flat plate area - Enter the aircraft equivalent flat plate area in square feet. 
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19. vertical projected area - Enter the vertical projected area of the aircraft in 
square feet. 

20. horizontal tail area - Enter the horizontal tail area in square feet. If no 
horizontal tail, enter a 0. 

21. horizontal tail span - Enter the horizontal tail span in feet. 

22. horizontal tail Cl - Enter the expected lift coefficient for the horizontal 
tail. 

23. horizontal tail Cdo - Enter the horizontal tail profile drag coefficient. 

24. vertical tail area - Enter the area of the vertical tail in square feet. If no 
vertical tail, enter a 0. 

25. vertical tail span - Enter the span of the vertical tail in feet. 

26. vertical tail Cl - Enter the expected lift coefficient for the vertical tail. 

27. vertical tail Cdo - Enter the vertical tail profile drag coefficient. 

28. 2P longitudinal input - Enter the 2/revolution fixed cyclic input. 

29. 3P input - Enter the 3/revolution fixed cyclic input. Enter 0 for no 3P 
input. Regardless of 3P input you must enter an offset value between 0 
and 119 degrees. 

30. auxiliary thrust - Enter the expected auxiliary thrust in pounds. 

31. wing area - Compound helicopter edit menu only. Enter the wing area in 
square feet. 

32. wing span - Compound helicopter edit menu only. Enter the wing span in 
feet. If you included the segment of fuselage encompassed by the wing in 
the wing area value include it in the span value also. 

33. wing Cl - Compound helicopter edit menu only. Enter the expected lift 
coefficient for the wing. 

34. wing Cdo - Compound helicopter edit menu only. Enter the wing profile 
drag coefficient. 
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35. wing efficiency factor - Compound helicopter edit menu only. Enter the 
wing’s expected efficiency factor. 

If you have chosen Compound Helicopter Analysis you will be prompted to 
answer three additional questions: 1. Do you want to fix the Tip Path Plane? 2. Do you 
want JANRAD to set thrust equal to total drag? 3. Do you want JANRAD 2P_RVRto 
use RVR based 2/rev Cyclic Pitch? It is recommended that the tip path plane be fixed for 
compound helicopters. If you chose to fix the angle you will be prompted to enter the 
angle in degrees. If you chose not fix the angle, the program will modulate the angle to 
achieve desired performance characteristics. If you chose not to set auxiliary thrust to 
total drag JANRAD use whatever value was entered in the edit menu. Only chose to use 
2/rev cyclic pitch if you are performing high speed (greater than 190 kts) slowed rotor 
analysis. If Conventional Helicopter Analysis was selected, these questions are bypassed 
and you are taken directly to the save instructions page. 

After all the required parameters are entered, you will be prompted to save the data 
under a user designated file name. The file name may be up to twenty characters long. 
JANRAD 2P_RVR will save the file to the current directory as the file name with a 
“.mat” extension. If you were editing a data file, press <ENTER> if you want to save the 
file with the original file name. Eigure 48 shows the Save Instructions page. 


*** SAVE IH3TRUCTI0H3 *** 

A. Save the new data to a specificed file name. 

B. Dg not use an extension or quotations. 

C. Use the letter/number combinations of 20 characters or less. 

D. The file will be saved with a ".mat" extension. 

ex: dsgn_2 

E. If you made no changes, press < Enter >, the file will 
be saved with the original name. 

save file as: | 

Eigure 48. Save Instructions Page. 
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The Execution Menu page is displayed next. Enter 1 if you want to proceed to 
next step in Rotor Performance Analysis, a 2 if you would like to go back and edit your 
data further or if you want to create a new file. To quit the program, select 3. 

After selecting Rotor Performance Analysis, JANRAD 2P_RVR will ask you 
several more questions based on user defined aircraft design inputs. Each question is 
discussed below. 

If a Compound Helicopter Analysis was selected, JANRAD 2P_RVR will display 
the percent of aircraft weight carried by the wing. You may choose to alter the 
percentage of weight carried by the wing or leave it at the original value. If you decide to 
change the value the program will calculate and display the new required wing area. 

Regardless of Conventional or Compound analysis JANRAD will check to 
deter mi ne if hover/low speed (airspeed less than 10 kts) performance analysis is about to 
be performed. If hover/low speed analysis has been selected you will be asked if there 
are any auxiliary vertical thrusters, i.e. direct jet thrusters or lift fans, installed on your 
aircraft. If installed, enter the amount of thrust in pounds provided by these devices. 

If a Compound Helicopter Analysis was selected and you chose not to set 
auxiliary thrust to total drag, JANRAD 2P_RVR will check to deter mi ne if a thrust deficit 
for you desired flight condition exist. If a deficit exists, it will be displayed and you will 
be given another chance to set auxiliary thrust equal to total drag. If you decide to 
continue with the original auxiliary thrust value, a warning message is displayed and the 
program proceeds to the next analytical step. 

If you previously chose to employ 2/rev cyclic pitch, JANRAD 2P_RVR will 
display the current value of the advance ratio, p, and if required, suggest altering p to the 
optimum value for RVR operation. Advance ratio must be greater than or equal to 1.0 to 
enable RVR based cyclic pitch. If you allow the program to change p the rotational 
velocity of your rotor will be changed to reflect the optimized p. 
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After all these inputs have been satisfied, the program will display an ‘ANALYSIS 
IN PROGRESS’ message and several status messages. If your inputs fail to produced a 
converged solution using either conventional or RVR trim criteria the message shown in 
Figure 49 will be displayed. 

This configuration will not trim 
Try a lower airspeed or a new design 

Program execution terminated 

??? *** EWD OF PROGRAM *** 

» _ 

Figure 49. Program Termination Page. 


From this page you must restart JANRAD 2P_RVR from the command prompt. 
A ‘ROTOR TRIMMED’ message will be displayed when the analysis is complete. 

2. Output 

Upon completion of the analysis, JANRAD 2P_RVR will ask if you want the 
performance results displayed on screen. If you answer ‘yes’ two screens displaying 
input data and two screens displaying output data will be displayed. The final output 
screen is displayed in Figure 50. 



CALCULATED DATA 

COHTINUED 


uh60a3 




solidity (sigma) 

0.082 



Disk Loading 

0.00 

lbs/ft'"2 


Figure of Merit = 

0.00 



CT/sigma 

0.081 



CQ/sigma 

0.0079 



CH/sigma 

0.0000 


Tip Mach numher of adv blade= 

0.808 



Advance ratio = 

0.285 


Rotor 

thrust reguired (TPP)= 

15366 

lbs 


Rotor power reguired= 

1984 

h.p. 


Rotor torgue= 

41243 

ft-lbs 

press 

any key to continue... 




Figure 50. Fourth Output Page. 
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Each of these screens is displayed one at a time. As soon as you advance to the 
next screen, you cannot return to the previous screen. A copy of the four screens is saved 
to the current directory with the file name you entered at the Save Instructions page and 
concatenated with a ‘.prf’ extension. If you answer ‘no’ the four screens will still be 
saved and you will be asked if you would like to view any of the JANRAD 2P_RVR 
created plots. If you answer ‘yes’ the Available Plots page, seen below in Figure 51, will 
be displayed. 


**** Available Plots **** 

1. Rotor Asimuth Angle v. Cyclic Pitch 

2 . Blade AEimuth Position v. Airload 

3. Rotor Airload Distribution 

4. Blade Azimuth Position v. Lift 

5. Rotor Lift Distribution 

6. Angle of Attack Distribution 

7. Rotor Tip Path Plane Stall pattern and Thrust Distribution 

8. Rotor Trim Lissajous Figure 

9. Lift Distribution (Contour Plot] 

10. CL Distribution (Contour Plot) 

11. Hach Humber Distribution (Contour Plot) 

Please choose a plot by entering 1-11.| 


Figure 51. Available Plots Page. 


Selected plots will be displayed in new windows labeled with the appropriate 
figure numbers. These plots may be simply viewed on the screen or printed by choosing 
the appropriate MATFAB menu bar commands. If they are viewed on screen it is 
recommended that you maximize the window for best clarity and plot resolution. When 
printing, set your printer to landscape mode for the best results. After you have 
completed viewing the selected figure, either close or minimize the window. A message 
will be in the command window asking if you want to view another plot. Answer ‘yes’ 
or ‘no’. If you chose ‘yes’ simply repeat the process until you have viewed all required 
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plots. If you choose to view plots 2, 4, or 7, you must manually position the plot titles. 
To do this, simply use the mouse to position the cursor at the point where you would like 
the title to begin and click the left button. Examples of each plot in Figures 52-62. 



Figure 52. Cyclic Pitch versus Rotor Azimuth Angle. 
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Figure 53. Blade Position versus Air Load. 



Fite Edit View Insert Tools Window Help 


Figure No. 3 




DcSHS \ ^ /< / ^ ^ O 

Airload Distribution At 120.1557 Ws, advance ratio=0.28189,2P inpul=4 degs, 3P input=0 degs, 3P offset=0 degs 


Figure 54. Air Load Distribution. 
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Figure 55. Blade Position versus Lift. 



File Edt View Insert Tools Window Help 


Figure No. 5 




iDcsgs K A y ^ ^ o 

Lifl Distribution At 120.1557 Kts, advance ratio=0.281B9,2P input=4 degs, 3P input=0 degs. 3P offset=0 degs 


Figure 56. Lift Distribution. 
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Figure 57. Angle of Attack Distribution. 



Figure 58. Actual Versus Ideal Rotor Tip Path Plane Angle and Thrust Distribution. 
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Figure 59. Lissajou Figure. 



Figure 60. Air Load Contour Plot. 
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Figure No. 10 


BBfEil 


RIe Edit View Insert Tools Window Help 

Dcsy# / ^ ^ o 



0 


Figure 61. Cl Distribution Contour Plot. 



Figure 62. Mach Number Distribution. 


57 































If you chose ‘no’, the Performance Output Data Instructions page will be 
displayed. See Figure 63. 

*** PERFCiPHAIICE OUTPUT DATA IHSTRUCTI0H3 *** 

A. Single value data saved to a file named: 

"ulieOaU.ptf" 

B. This is a text file, use the tyoe command to view the 
file OE use a text editor to view/print the file. 

A duplicate file name exists, or the file 
cannot be found. 

C. HatEix and vector data saved to a file named: 

"uh60a3_p.mat" 

D. This is a ".mat" binary file, use the load command 
to retrieve the data for plotting. 

*** EHD OF PERFORHAUCE *** 

press any key to continue... 

Figure 63. Performance Output Data Instructions Page. 

3. Saved Data 

The vector and matrix data generated by the rotor performance module of 
JANRAD 2P_RVR is saved to the current directory with the file name you entered at the 
Save Instruction page and concatenated with a ‘_p’ and a ‘.mat’ extension. The 
following vectors and matrices are saved. 

1. r - a vector containing the radii (in feet) to the center of each blade 
element. This is a row vector containing n elements, where n equals the number of blade 
elements plus one. 

2. ^ - a vector containing the azimuth angles (in degrees) around the rotor 
disk. This is a column vector containing n elements, where n equals the number of 
azimuth sectors. 

3. V, - a vector containing the induced velocity distribution (in ft/sec)along 
the rotor blade. This is a row vector containing n elements, where n equals the number of 
blade elements plus one. 
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4. 6* - a vector containing the 1/rev cyclic pitch (in degrees) with respect to 
azimuth angle at 0.7 r/R. This is a column vector containing n entries, where n equals the 
number of azimuth sectors. 

5. ^xp ~ ^ vector containing the total cyclic pitch, 1/rev + 2/rev +3/rev, (in 

degrees) with respect to azimuth angle at 0.7 r/R. This is a column vector containing n 
entries, where n equals the number of azimuth sectors. 

6. P, - ^ vector containing blade pitch (in degrees) along the rotor blade. 

This is a row vector containing n elements, where n equals the number of blade elements 
plus one. 

7. a - a matrix containing the angle of attack (in degrees) distribution. The 
matrix is with m rows and n columns, where m equals the number of azimuth sectors and 
n equals the number of blade elements plus one. 

8. ~ ^ vector containing the total thrust (in pounds) at each azimuth sector. 

This is a column vector containing n elements, where n equals the number of azimuth 
sectors. 

9. - a vector containing the total thrust moment (in ft-lbs) at each 

azimuth sector. This is a column vector containing n elements, where n equals the 
number of azimuth sectors. 

10. DM^ - a vector containing the total drag moment (in ft-lbs) at each 

azimuth sector. This is a column vector containing n elements, where n equals the 
number of azimuth sectors. 

11. dT- a matrix containing the differential thrust (in pounds) distribution. 
This is a matrix with m rows and n columns, where m equals the number of azimuth 
sectors and n equals the number of blade elements plus one. 

12. dM - a matrix containing the differential thrust moment (in ft-lbs) 
distribution. This is a matrix with m rows and n columns, where m equals the number of 
azimuth sectors and n equals the number of blade elements plus one. 
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13. dD - a matrix containing the differential drag (in pounds) distribution. 
This is a matrix with m rows and n columns, where m equals the number of azimuth 
sectors and n equals the number of blade elements plus one. 

14. Ciplot - a matrix containing the CL distribution over the rotor disk. This 
is a matrix with n rows and n columns, where n equals the number of blade elements. 

15. Coplot - a matrix containing the CD distribution over the rotor disk. This 
is a matrix with n rows and n columns, where n equals the number of blade elements. 

When JANRAD 2P_RVR has completed saving the data, the Execution Menu 
page will be displayed. See Figure 64. Choose the desired function. 

*** EXECUTION MENU *** 

1. Rotor Performance Analysis 

2. Change Data 

3. Quit 

Enter a 1, 2, or 3 »| 

Figure 64. Execution Menu Page. 

If you chose option 2, you will be able to reload the file you were just working on 
or load an entirely different file. If you chose to load an entirely different file, it must be 
of the same type (either Conventional or Compound Helicopter) as the one you just 
completed working on. If you need to change the type of helicopter (conventional or 
compound) you must select option 3 from the Execution Menu page and restart JANRAD 
2P_RVR. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

The original intent of this investigation was to deter mi ne to what extent helicopter 
forward flight could be improved using Reverse Velocity Rotor (RVR) technology and 
Higher Harmonic Stall Control (HHSC). JANRAD 2P_RVR computational analysis 
shows that HHSC alone can produce increases in forward flight speeds of up to ten 
percent in conventional helicopters. Although these improvements were not as radical as 
Arcidiacono predicted they are still significant. The redistribution of the lift over the 
rotor disk allows the helicopter to operate with a greater percentage of the retreating side 
of the rotor system engulfed in reverse flow while still meeting the rigid trim criteria 
imposed by the program’s original trim module. Employment of the 2/rev feathering 
incurred minimal power penalties and allowed the blade thrust coefficient to approach 
ideal values. Use of the HHSC worked equally well when used with compound 
helicopter configurations up to 175 knots, although the extra weight and drag caused by 
the wing increased power requirements. 

At typical helicopter airspeeds (less than 160 knots) the RVR is somewhat 
inefficient. The RVR’s basic shape does not perform nearly as well as the exotic blades 
currently in use by rotorcraft but, the use of a symmetrical, double ended airfoil allows 
the generation of lift with both forward and reverse flows across the blades, which would 
allow the generation of lift on the retreating side at a dramatically wider range of 
airspeeds than conventionally shaped blades. RVR equipped rotorcraft would require 
more power to hover and climb than a rotorcraft equipped with proprietary non-linear 
twist or taper. The true potential of the RVR is best exploited at high forward flight 
speeds. As flight speed increases the rotor must be slowed based on predetermined 
optimum advance ratio scheduling. As demonstrated in the brief JVHL study, the RVR is 
most efficiently placed in this environment when it is used on an aircraft equipped with a 
compound wing and auxiliary thrust devices. An RVR equipped aircraft would spend the 
majority of its time in cruise flight, minimizing time spent hovering and at low airspeeds. 
At cruise speed the RVR must be controlled by a unique RVR specific 2/rev HHSC 

control system. This system builds upon the conventional 2/rev HHSC system and is 
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required to control the angle of attack of the blades around the rotor system’s azimuth to 
maxi miz e the RVR’s unique performance. This angle of attack manipulation allows for 
lift generation in all quadrants of the system at cruise speed. The predictable lift 
generation allows for controlled aerodynamic flight of the rotor disk. 

B. RECOMMENDATIONS 

Improvements to the program can be made in several ways. First, the graphical 
user interface features within MATLAB should be harnessed to allow users to enter data 
via editable dialogue and text boxes. The program can also be modified to allow for 
offsetting the center of gravity from the center of the rotor hub during the trimming 
process. As more data are made available more airfoils should be added to the selection 
list to allow for broader experimentation. 

Although JANRAD 2P_RVR has the capacity to accommodate 3/rev inputs, the 
study was limited to 2/rev inputs. Investigation into the effects of combining 3/rev inputs 
with both a non-RVR and RVR system may yield unforeseen benefits or behaviors. 

The development of a HHSC/RVR blade dynamics and stability and control 
module should be considered. The module should include provisions for both analytical 
and visual output. These modules should also be able to ‘plug-in’ to the existing 
JANRAD 2P_RVR performance module to give the user a complete preliminary design 
tool. 

With the advent of new V/STOL technologies the use of lift fans, direct jet 
thrusters, and other lift augmentation devices will become more commonplace. Future 
versions of the program must be able to fully incorporate these devices and should be 
able to integrate their effects with the traditional helicopter lift (main rotor) and anti¬ 
torque (tail rotor) devices. New techniques will have to be developed to combine the 
effects of these systems to accurately determine downwash velocities, climb 
performance, and controllability. 
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APPENDIX A. JANRAD.M 


% Joint Army Navy Rotorcraft Analysis and Design 
% JANRAD 
% Version 2P_RVR 

% All modification to original code written by 
% MAJ Steven Van Riper 
% August 2003 
% 

% Core programming technology 
% developed by 
% E. R. Wood, Ph.D. 

% MAJ Bob Nicholoson 
% MAJ Walter Wirth 
% September 1993 
% 

% This program is an interactive preliminary design tool 
% developed to aid the design student in determination of 
% initial rotorcraft configurations and in the calculation 
% of performance and other parameters with or without 
% IP + 2P + 3P rotor control. 

% The program will work for conventional or compound rotorcraft. 
% The program can take into account various types of powered lift 
% including lift fans and direct jet thrusters. 

% It will provide accurate data for airspeeds less than 10 
% knots and greater than or equal to 50 knots. 

% Variable list for janrad.m, trim.m, thrcalc.m, tmcalc.m, dmcalc.m, 
% hh02clcd.m, vrl2clcd.m, plotter.m, and perf.m 

% a lift curve slope of rotor system airfoil 

% Adisk area of rotor disk 
% Afh fuselage equivalent flat plate drag area 
% Afv vertical projected area (fuselage area under disk) 

% airfoil rotor system airfoil type (HH02WR12) 

% alpha angle of attack, rotor blade radial segment 
% alpham matrix of alpha values for each blade element 
% alphaT rotor tip path plane angle 

% alphaTfix determines if user fixes alphaT to specified value 
% area area of Lissajous figures in plotter.m 
% b number of rotor blades 

% B tip loss parameter 

% betao rotor coning angle 

% betat geometric angle, rotor blade radial segment 

% bhoriz span, horizontal tail 

% bvert span, vertical tail 

% bwing span, wing 
% cblade chord, rotor blade 

% CD drag coefficient, rotor blade radial segment 

% CDohoriz profile drag coefficient, horizontal tail 
% CDovert profile drag coefficient, vertical tail 
% CDOwing profile drag coefficient, wing 
% CDhoriz drag coefficient, horizontal tail 
% CDplot matrix of CD values for each blade element 
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% CDvert drag coefficient, vertical tail 
% CDwing drag coefficient, wing 
% CH rotor H force coefficient 

% CH_sig CH/solidity 

% CL lift coefficient, rotor blade radial segment 

% CLhoriz lift coefficient, horizontal tail 
% CLplot matrix of CL values for each blade element 
% CLvert lift coefficient, vertical tail 
% CLwing lift coefficient, wing 
% CQ rotor torque coefficient 

% CQ_sig CQ/solidity 

% CT rotor thrust coeeficient 

% CT_sig CT/solidity 

% dD differential drag, rotor blade radial segment 

% ddD differential drag, rotor blade tip 
% ddDM differential drag moment, rotor blade tip 

% ddM differential thrust moment, rotor blade tip 
% ddT differential thrust, rotor blade tip 
% delM change in total thrust moment 
% dev AC deviation used within Adjusting Collective routine 

% devCY deviation used within Adjusting Cyclic routine 

% devMT deviation used within Mean Thrust routine 

% devTC deviation used within Trimming Collective routine 
% Dftotal resultant of fuselage drag and aux thrust 
% Dfuse total drag generated by no rotor bodies 
% DL disk loading 

% dM differentail thrust moment, rotor blade radial segment 

% DMpsi total blade drag moment at specific azimuth angle 
% dr rotor blade radial segment width 

% Drotor rotor system drag 

% dT differential thrust, rotor blade radial segment 

% Dtplot matrix containing dTs for all blade elements 

% Dhoriz drag, horizontal tail 

% dthetadM change in cyclic pithch with change in thrust moment 
% Overt drag, vertical tail 

% Owing drag, wing 
% e effective hinge offset 

% ewing wing efficiency factor 

% filename name of input file 
% FM figure of merit 

% grip lenght of inner non-aerodynamic protion of blade 
% GW aircraft gross wieght 

% helochoice determines choice of compound or conventional helo 
% Hrotor rotor H force 

% hub gives size of hub in plotter.m 

% lamdaT forward flight induced velocity parameter 
% Lftotal total lift generated by non-rotor bodies 
% Lhoriz lift, horizontal tail 

% lifterror difference between actual lift and desired lift 
% LoverT Lift over thrust percentage 

% Lvert lift, vertical tail 

% Lwing lift, wing 

% Mlc first harmonic (cosine) thrust moment coefficient 
% Mis first harmonic (sine) thrust moment coefficient 
% Machtip Mach number at rotor blade tip 

% mblade mass of rotor blade 
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% Mpsi total blade thrust moment at specific azimuth angle 
% mu advance ratio 

% muchoice determines if user would like JANRAD to change mu to 
% a desired value 

% muerror difference between actual mu and desired mu 
% naz number of azimuth sectors 
% nbe number od blade elements 

% omega rotor rotational velocity 

% PA pressure altitude 

% phi inflow angle, rotor blade radial segment 

% phitip inflow angle, rotor blade tip 

% plotchoice determine which plot to display in plotter.m 

% Protor power required by the rotor 

% psi azimuth angle 

% q dynamic pressure 

% Qrotor rotor torque 

% r radius, rotor blade radial segment 

% R rotor blade radius 

% Rbar Reff-e 

% RbarT rT*Rbar 

% Reff effective rotor radius, tip loss 

% rho ambient air density 

% rT location of resultant thrust vector 

% rvrCP determines if normal or rvr cyclic pitch scheduling is used 

% solidity solidity 

% Shoriz area, horizontal tail 

% Svert area, vertical tail 

% Swing area, wing 

% T rotor thrust 

% Taux auxiliary thrust 

% Tauxchoice determines if user would like JANRAD to match Taux 

% requirements 

% temp ambient air temperature 

% theta cyclic pitch 

% thetahub hub thetas in plotter.m 

% theta Ic first harmonic (cosine) of cyclic pitch 

% thetals first harmonic (sine) of cyclic pitch 

% theta2c second harmonic (cosine) of cyclic pitch 

% thetaSs third harmonic (sin) of cyclic pitch 

% thetaJPoffset 3/rev offset 

% thetao collective pitch at .7 r/R 

% thetaXP total vine of cyclic pitch with all inputs 

% Tpsi total blade thrust at specific azimuth angle 

% twist geometric rotor blade twist 

% Up vertical component of velocity 

% Uptip vertical component of velocity at tip 

% Ut horizontal component of velocity 

% Uttip horizontal component of velocity at tip 

% vertaux vertical auxiliary thrust in pounds 

% vertauxchoice determines if auxiliary vertical thrust devices are 

% installed on the aircraft 

% vi induced velocity 

% Vinf forward airspeed 

% Vtip tip speed 

% wblade wieght of rotor blade 

% wingchoice determines if user wnats to allow JANRAD to resize the wing 
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% wpercent percent of aircraft gross wieght to be carried by the wing 

% clearing all the variables is the MATLAB environment 

clear all 

clc 

disp ('') 
disp ('') 
disp 


disp C 

* 



disp C 

* 

Joint Army/Navy Rotorcraft Analysis and Design 

disp C 

* 

(JANRAD) 


disp C 

* 



disp C 

* 

Version 2P_RVR 


disp C 

* 

October 2003 


disp C 

* 



disp 





disp ('') 

% 

% Determine if the user wants to analyze conventional or compound 
% helicopter models 

% 

flag=l; 
while flag > 0 
dispC) 

dispCDo you want to analyze conventional or compound rotorcraft? ') 
dispC) 

helochoice=input(' 1. conventional or 2. compound »'); 

while isempty(helochoice), 
dispC ') 

dispC You must enter a 1 or 2') 

dispCDo you want to analyze conventional or compound rotorcraft? ') 
helochoice=inputCl. conventional or 2. compound »'); 
end 

if helochoice~= 1 & helochoice~=2 
dispC ') 

dispC *** Enter a 1 or 2 ***') 
else 
flag=0; 
end 

end 

clc 

checkl=l; 
while check 1 > 0 
checkl=l; 
dispC ') 

dispCDo you want to edit any existing file or create a new one?') 

dispC') 

check=1; 

while check > 0 

answer=inputC 1. edit existing file 2. create new file »'); 


% *** If editing an existing file: get file name, display edit 



% menu, allow changes to selected variables, and save under 
% desired file name. Loads to and saves from current 
% directory a a .mat file. *** 

if answer== 1 
clc 

dispC ') 
dispC ') 

dispC * * * LOAD INSTRUCTIONS * * *') 

dispC ') 

dispC A. Input the name of the file to edit.') 

dispC B- The file was saved in your previous session') 

dispC with a ".mat" extension.') 

dispC C. Do not include the extension or quaotations.') 

dispC ') 

dispC ex: dsgnl') 

flag=0; 

while flag < 1 

filename l=inputC name of input file: ','s'); 
e val( ['flag=existC",filename 1,'.mat");']) 
if flag < 1 
disp ('') 

disp (' The file does not exist, try again or <Ctrl-C>') 
disp (' to exit program.') 
end 
end 

eval(['load ',filenamel]) 
while check > 0 
if helochoice==l 
clc 

dispC ') 

dispC * * * CONVENTIONAL HELICOPTER EDIT MENU * * *') 

dispC ') 

dispC 1 • pressure altitude 2. temperature') 

dispC 3. airspeed 4. gross weight') 

dispC 5. number of blades 6. blade radius') 

dispC 7. blade chord 8. hinge offset') 

dispC 9. blade grip length 10. blade twist') 

dispC 11. blade wieght 12. number of blade elements') 

dispC 13. rotational velocity 14. # azimuth sectors') 
dispC 15. lift curve slope 16. airfoil') 

dispC 17. collective pitch 18. flatplate area') 

dispC 19. vert projected area 20. horizontal tail area') 

dispC 21. horizontal tail span 22. horizontal tail CL) 

dispC 23. horizontal tail CDo 24. vertical tail area') 
dispC 25. vertical tail span 26. vertical tail CL') 

dispC 27. vertical tail CDo 28. 2P input (cosine)') 

dispC 29.3P input ') 

dispC 0. NO CHANGES') 

alphaTfix=0;itercount=0;helochoice=l;rvrCP=2; 
else 
clc 

dispC ') 

dispC *** COMPOUND HELICOPTER EDIT MENU ***') 
dispC ') 

dispC 1 • pressure altitude 2. temperature') 
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dispC 

3. airspeed 4 

gross weight') 

dispC 

5. number of blades 

6. blade radius') 

dispC 

7. blade chord 

8. hinge offset') 

dispC 

9. blade grip length 

10. blade twist') 

dispC 

11. blade wieght 

12. number of blade elements') 

dispC 

13. rotational velocity 

14. # azimuth sectors') 

dispC 

15. lift curve slope 

16. airfoil') 

dispC 

17. collective pitch 

18. flatplate area') 

dispC 

19. vert projected area 

20. horizontal tail area') 

dispC 

21. horizontal tail span 

22. horizontal tail CL') 

dispC 

23. horizontal tail CDo 

24. vertical tail area') 

dispC 

25. vertical tail span 

26. vertical tail CL') 

dispC 

27. vertical tail CDo 

28. 2P input') 

dispC 

29. 3P input 30. auxiliary thrust') 

dispC 

31. wing area 

32. wing span') 

dispC 

33. wing CL 

34. wing CDo') 

dispC 

35. wing efficiency factor') 

dispC 

0. NO CHANGES') 



alphaTfix=0;itercount=0;helochoice=2;rvrCP=2; 

end 

choice=input(' Select the parameter to change; '); 
if isempty(choice), 
choice=0; 
end 

if choice==l, 
clc 

dispC ') 

PA 

tmp=PA; 

PA=input('Pressure altitude (ft);'); 
if isempty(PA), 

PA=tmp; 

end 

elseif choice==2, 
clc 

dispC ') 
temp 

tmp=temp; 

temp=input('temperature (deg F):'); 
if isempty(temp)„ 
temp=tmp; 
end 

elseif choice==3, 
clc 

dispC ') 

Vinf=Vinf/1.69 

tmp=Vinf; 

Vinf=input('Airspeed (knots): ')*1.69; 
if isempty(Vinf), 

Vinf=tmp*1.69; 

end 

elseif choice==4, 
clc 

dispC ') 

GW 

tmp=GW; 
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GW=input('Aircraft Groos weight (lbs): '); 
if isempty (GW), 

GW=tmp; 

end 

elseif choice==5, 
clc 

dispC ') 
b 

tmp=b; 

b=input('Number of blade: '); 
if isempty(b), 
b=tmp; 
end 

elseif choice==6, 
clc 

dispC ') 

R 

tmp=R; 

R=input('Blade radius; center of hub to blade tip (ft): '); 
if isempty(R), 

R=tmp; 

end 

elseif choice==7, 
clc 

dispC ') 

cblade 

tmp=cblade; 

cblade=input('Blade chord (ft);'); 
if isempty(cblade), 
cblade=tmp; 
end 

elseif choice==8, 
clc 

dispC ') 
e 

tmp=e; 

e=inputCHinge offset (ft);'); 
if isempty(e), 
e=tmp; 
end 

elseif choice==9, 
clc 

dispC ') 
grip 

tmp=grip; 

grip=inputCnon-aerodyn inboard portion of blade (ft): '); 
if isempty (grip), 
grip=tmp; 
end 

if grip < le-1, 
grip=le-10; 
end 

elseif choice==10 
clc 

dispC ') 

twist=-twist*57.3 
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tmp=twist; 

twist=input('blade twist (deg): '); 
if isempty (twist), 
twist=abs(tmp)/57.3; 
else 

twist=abs(twist)/57.3; 

end 

elseif choice==l 1, 
clc 

dispC ') 
wblade 
tmp=wblade; 

wblade=input('Weight of aero portion of blade (lbs): '); 
if isempty(wblade), 
wblade=tmp; 
end 

elseif choice==12, 
clc 

dispC ') 
nbe 

tmp=nbe; 

nbe=input('number of blade elements: '); 
if isempty(nbe), 
nbe=tmp; 
end 

elseif choice==13, 
clc 

dispC ') 
omega 
tmp=omega; 

omega=input('Rotor rotational velocity (rad/sec): '); 
if isempty (omega), 
omega=tmp; 
end 

elseif choice==14, 
clc 

dispC ') 
naz 

tmp=naz; 

naz=inputCnumber of azimuth sectors: '); 
if isempty (naz), 
naz=tmp; 
end 

elseif choice==15, 
clc 

dispC ') 
a 

tmp=a; 

a=inputCAverage lift curve slope (CL vs alpha);'); 
if isempty(a), 
a=tmp; 
end 

elseif choice==16, 
clc 

dispC') 

airfoil 
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tmp=airfoil; 
flag=l; 
while flag > 0 

airfoil=input('Airfoil 1. HH-02 2. VR-12 3. SC1095 4. 0012 »'); 
if isempty( airfoil), 
airfoil=tmp; 
end 

if airfoil==l, 
flag=0; 

elseif airfoil==2, 
flag=0; 

elseif airfoil==3, 
flag=0; 

elseif airfoil==4, 
flag=0; 
else 
dispC ') 

dispC *** Enter a 1, 2, 3, or 4 *** ') 
end 
end 

elseif choice==16 & helochoice==2, 
clc 

dispC") 
airfoil 
tmp=airfoil; 
flag=l; 
while flag > 0 

airfoil=input('Airfoil 1. HH-02 2. VR-12 3. SC1095 4. 0012 5. RVR»'); 
if isempty(airfoil), 
airfoil=tmp; 
end 

if airfoil==l, 
flag=0; 

elseif airfoil==2, 
flag=0; 

elseif airfoil==3, 
flag=0; 

elseif airfoil==4, 
flag=0; 

elseif airfoil==5, 
flag=0; 
else 
dispC ') 

dispC *** Enter a 1, 2, 3, 4, or 5 *** ') 
end 
end 

elseif choice==17, 
clc 

dispC ') 

thetao=thetao*57.3 

tmp=thetao; 

thetao=inputCcollective pitch at .7 r/R (deg): ')/57.3; 
if isempty(thetao), 
thetao=tmp/57.3; 
end 

elseif choice==18, 
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clc 

dispC ') 

Afh 

tmp=Afh; 

Afh=input('Aircraft equivalent flatplate area (ft^2):'); 
if isempty(Afh), 

Afh=tmp; 

end 

elseif choice==19, 
clc 

dispC ') 

Afv 

tmp=Afv; 

Afv=input('Verical projected area (ft^2):'); 
if isempty(Afv), 

Afv=tmp; 

end 

elseif choice==20, 
clc 

dispC ') 

Shoriz 

tmp=Shoriz; 

Shoriz=input('Horizontal tail area (ft''2): '); 
if isempty(Shoriz), 

Shoriz=tmp; 

end 

if Shoriz < le-10, 

Shoriz=le-10; 

end 

elseif choice==21, 
clc 

dispC ') 

bhoriz 

tmp=bhoriz; 

bhoriz=inputCHorizontal tail span (ft): '); 
if isempty(bhoriz), 
bhoriz=tmp; 
end 

if bhoriz < le-10, 
bhoriz= le-10; 
end 

elseif choice==22, 
clc 

dispC ') 

CLhoriz 

tmp=CLhoriz; 

CLhoriz=inputCExpected CL for the horizontal tail: '); 
if isempty (CLhoriz), 

CLhoriz=tmp; 

end 

elseif choice==23, 
clc 

dispC ') 

CDohoriz 

tmp=CDohoriz; 

CDohoriz=inputCHorizontal tail prolfile drag coef (CDo): '); 

72 



if isempty(CDohoriz), 

CDOhoriz=tmp; 

end 

elseif choice==24, 
clc 

dispC ') 

Svert 

tmp=Svert; 

Svert=input('Vertical tail area (ft^2): '); 
if isempty(Svert), 

Svert=tmp; 

end 

if Svert < le-10, 

Svert=le-10; 

end 

elseif choice==25, 
clc 

dispC ') 
bvert 

tmp=bvert; 

bvert=input('Vertical tail span (ft): '); 
if isempty(bvert), 
bvert=tmp; 
end 

if bvert < le-10, 
bvert= le-10; 
end 

elseif choice==26 
clc 

dispC ') 

CLvert 

tmp=CLvert; 

CLvert=input('Expected CL for the vertical tail: '); 
if isempty(CLvert), 

CLvert=tmp; 

end 

elseif choice==27, 
clc 

dispC ') 

CDovert 

tmp=CDovert; 

CDovert=inputCVertical tail profile drag coef (CDo);'); 
if isempty(CDovert), 

CDovert=tmp; 

end 

elseif choice==28, 
clc 

dispC ') 

twoPinput 

tmp=twoPinput; 

twoPinput=inputCLongitudinal 2P harmonic input (degs): '); 
if isempty (twoPinput), 
twoPinput=tmp; 
end 

elseif choice==29 
clc 
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disp(") 

threePinput 

tmp=threePinput; 

threePinput=input('3P harmonic input (degs): '); 

dispC) 

threePoffset 

tmp=threePoffset; 

threePoffset=input('3P input offset, must be between 0 and 119 degs 
if isempty (threePinput) 
threePinput=tmp; 
end 

elseif helochoice==2 & choice==30, 
clc 

dispC ') 

Taux 

tmp=Taux; 

Taux=input('Auxiliary thrust (lbs): '); 
if isempty (Taux), 

Taux=tmp; 

end 

elseif helochoice==2 & choice==31, 
clc 

dispC ') 

Swing 

tmp=Swing; 

Swing=input('Wing Area (ft''2):'); 
if isempty (Swing), 

Swing=tmp; 

end 

if Swing < le-10 
Swing=le-10; 
end 

elseif helochoice==2 & choice==32, 
clc 

dispC ') 

bwing 

tmp=bwing; 

bwing=inputCWing span (ft);'); 
if isempty(bwing), 
bwing=tmp; 
end 

if bwing < le-10, 
bwing= le-10; 
end 

elseif helochoice==2 & choice==33, 
clc 

dispC ') 

CLwing 

tmp=CLwing; 

CLwing=inputCexpected CL for the wing:'); 
if isempty(CLwing), 

CLwing=tmp; 

end 

elseif helochoice==2 & choice==34, 
clc 

dispC ') 


74 



CDowing 

tmp=CDowing; 

CDowing=input('Wing profile drag coef (CDo):'); 
if isempty(CDowing), 

CDowing=tmp; 

end 

elseif helochoice==2 & choice==35, 
clc 

dispC ') 

ewing 

tmp=ewing; 

ewing=input('Wing efficiency factor (e): '); 
if isempty(ewing), 
ewing=tmp; 
end 

if ewing < le-10, 
ewing=le-10; 
end 

elseif helochoice==2 & choice==0 
GW1=GW; 
clc 

ifVinf> 16.9 
disp(") 

dispCDo you want to fix the Tip Path Plane (recommended for compound 

helicopters),') 

alphaTfix=input('l. yes or 2. no?'); 
if alphaTfix==l 
dispC) 

alphaTinput=input('Enter the Tip Path Plane angle in degrees: '); 
alphaTinput=alphaTinput/57.3; 

dispC*** You have fixed the Tip Path Plane angle, please set auxiliary ***') 
dispC*** thrust to equal total drag to avoid erroneous solutions. ***') 
else 
dispC) 

dispCTip path plane will be determined by JANRAD') 
dispCPress any key to continue...') 
pause 
end 

% added 14-16 Sep 2003 to allow user to disable the retrim routine, select Taux based 
% on drag, and allow for RVR based cyclic inputs in trim.m 
flag=l; 
while flag > 0 
dispC) 

dispCDo you want JANRAD to set auxiliary thrust to equal total drag,') 
Tauxchoice=inputCl. yes or 2. no?'); 
while isempty(Tauxchoice), 
dispC') 

dispCYou must enter a 1 or 2') 

dispCDo you want JANRAD to set auxiliary thrust to equal total drag,') 
Tauxchoice=inputCl. yes or 2. no?'); 
end 

if Tauxchoice~=l & Tauxchoice~=2 
dispC') 

dispC *** Enter a 1 or 2 ***') 
elseif Tauxchoice==l 

dispCJANRAD will adjust auxiliary thrust to equal total drag') 
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dispCPress any key to continue...') 

pause 

flag=0; 

elseif Tauxchoice==2 

dispC*** Auxiliary thrust defined by user. ***') 
dispCPress any key to continue...') 
pause 
flag=0; 
else 
flag=0; 
end 
end 

rvrCP=2; 
if Vinf/1.69>200 
rvrCP=0; 
flag=l; 
while flag > 0 
disp(") 

disp('Do you want JANRAD to use RVR based (constrained) 2/rev Cyclic 

Pitch,') 

rvrCP=input('l. yes or 2. no?'); 
while isempty(rvrCP), 
disp('') 

disp('You must enter a 1 or 2') 

disp('Do you want JANRAD to use RVR based (constrained) 2/rev Cyclic 

Pitch,') 

rvrCP=input('l. yes or 2. no?'); 
end 

if rvrCP~=l & rvrCP~=2 
disp('') 

disp(' *** Enter a 1 or 2 ***') 
elseif rvrCP==l 
disp(") 

disp('RVR 2/rev Cyclic Pitch will be used') 
disp('Press any key to continue...') 
pause 
flag=0; 

elseif rvrCP==2 
disp(") 

disp('Conventional 2/rev cyclic pitch will be used') 
disp('Press any key to continue...') 
pause 
flag=0; 
else 
flag=0; 
end 
end 
else 
end 
else 
end 
clc 

dispC ') 
dispC ') 

disp(' *** SAVE INSTRUCTIONS ***') 

dispC ') 


76 



dispC A. Save the new data to a specificed file name.') 

dispC B. Do not use an extension or quotations.') 

dispC C. Use the letter/number combinations of 20 characters or less.') 

dispC D. The file will be saved with a ".mat" extension.') 

dispC ') 

dispC ex: dsgn_2') 
dispC ') 

dispC E. If you made no changes, press < Enter >, the file wilf) 
dispC be saved with the original name.') 

flag=l; 
while flag > 0 

filenameO=filename 1; 
filenamel=inputC save file as: ','s'); 
if isempty(filenamel) 
filename 1 =filenameO; 
end 

clear filenameO 
if length(filenamel) > 20, 
dispC ') 

dispC use 20 characters or less') 
flag=l; 
else 
flag=0; 
end 
end 

eval(['save ',filenamel]) 
check=0; 
elseif choice==0 
if helochoice== 1 
Swing=le-10; 
bwing=le-10; 

CLwing=0; 

CDowing=0; 

ewing=le-10; 

Taux=0; 

GW1=GW; 

else 

end 

clc 

dispC ') 
dispC ') 

dispC *** SAVE INSTRUCTIONS ***') 

dispC ') 

dispC A. Save the new data to a specificed file name.') 

dispC B- Do not use an extension or quotations.') 

dispC C. Use the letter/number combinations of 20 characters or less.') 

dispC D. The file will be saved with a ".mat" extension.') 

dispC ') 

dispC ex: dsgn_2') 
dispC ') 

dispC E. If you made no changes, press < Enter >, the file will') 
dispC be saved with the original name.') 

flag=l; 
while flag > 0 

filenameO=filename 1; 

filenamel=inputC save file as: ','s'); 
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if isempty(filenamel) 
filename 1 =filenameO; 
end 

clear filenameO 
if length(filenamel) > 20, 
dispC ') 

dispC use 20 characters or less') 
flag=l; 
else 
flag=0; 
end 
end 

eval(['save ',filenamel]) 
check=0; 

else 
dispC ') 

dispC Enter a displayed number ...press any key to continue') 
pause 
end 
end 


% *** If creating a new file; get input for required variables 
% and save under desired file name. Save to current 
% directory as a .mat file. *** 

elseif answer==2, 

alphaTfix=0;itercount=0; 

clc 

PA=inputCPressure Altitude (ft): '); 
while isempty(PA), 
dispC ') 

dispCYou must enter a numerical value') 
PA=inputCPressure altitude (ft): '); 
end 

temp=input('Temperature (deg F):'); 
while isempty(temp), 
dispC ') 

dispCYou must enter a numerical value') 
temp=input('Temperature (deg F): '); 
end 

Vinf=input('Airspeed (knots): ')*1.69; 
while isempty(Vinf), 
dispC ') 

dispCYou must enter a numerical value') 
Vinf=input('Airspeed (knots): ')*1.69; 
end 

GW=input('Aircraft gross wieght (lbs): '); 
while isempty(GW), 
dispC ') 

dispCYou must enter a numerical value') 
GW=input('Aircraft gross wieght (lbs):'); 
end 

b=input('Number of Blades: '); 
while isempty(b), 
dispC ') 
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dispCYou must enter a numerical value') 
b=input('Number of Blades: '); 
end 

R=input('Blade radius: center of hub to blade tip (ft): '); 
while isempty(R), 
dispC ') 

dispCYou must enter a numerical value') 

R=input('Blade radius: center of hub to blade tip (ft): '); 
end 

cblade=input('Blade chord (ft): '); 
while isempty(cblade), 
dispC ') 

dispCYou must enter a numerical value') 
cblade=input('Blade chord (ft): '); 
end 

e=input('Hinge offset (ft)'); 
while isempty(e) 
dispC ') 

dispCYou must enter a numercial value') 
e=input('Hinge offset (ft)'); 
end 

grip=input('Non-aerodynamic inboard portion of blade (ft): '); 
while isempty(grip) 
dispC ') 

dispCYou must enter a numercial value') 
grip=input('Non-aerodynamic inboard portion of blade (ft): '); 
end 

while grip < le-10, 
grip=le-10; 
end 

twist=input('Blade twist (deg): ')/57.3;,twist=abs(twist); 
while isempty(twist), 
dispC ') 

dispCYou must enter a numerical value') 
twist=input('Blade twist (deg): ')/57.3;,twist=abs(twist); 
end 

wblade=input('weight of aero portion of one blade (lbs): '); 
while isempty(wblade), 
dispC ') 

dispCYou must enter a numerical value') 
wblade=input('weight of aero portion if one blade (lbs): '); 
end 

nbe=input('Number of blade elements: '); 
while isempty(nbe), 
dispC ') 

dispCYou must enter a numerical value') 
nbe=input('Number of blade elements: '); 
end 

omega=input('Rotor rotational velocity (rad/sec):'); 
while isempty(omega), 
dispC ') 

dispCYou must enter a numerical value') 
omega=input('Rotor rotational velocity (rad/sec): '); 
end 

naz=input('Number of azimuth sectors: '); 
while isempty(naz), 
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dispC ') 

dispCYou must enter a numerical value') 
naz=input('Number of azimuth sectors: '); 
end 

a=input('Lift curve slope of rotor airfoil (CL vs alpha): '); 
while isempty(a), 
dispC ') 

dispCYou must enter a numerical value') 
a=input('Lift curve slope of rotor airfoil (CL vs alpha): '); 
end 

flag=l; 
while flag > 0 
if helochoice==l 

airfoil=input('Airfoil 1. HH-02 2. VR-12 3. SC1095 4. 0012 »'); 
while isempty(airfoil), 
dispC ') 

dispCYou must enter a numerical value') 
airfoil=input('Airfoil 1. HH-02 2. VR-12 3. SC1095 4. 0012 »'); 
end 

if airfoil~=l & airfoil~=2 & airfoil~=3 & airfoil~=4 
dispC ') 

dispC *** Enter a 1, 2, 3, or 4 ***') 
else 
flag=0; 
end 
else 

airfoil=input('Airfoil 1. HH-02 2. VR-12 3. SC1095 4. 0012 5. RVR »'); 
while isempty(airfoil), 
dispC ') 

dispCYou must enter a numerical value') 

airfoil=input('Airfoil 1. HH-02 2. VR-12 3. SC1095 4. 0012 5. RVR »'); 
end 

if airfoil~=l & airfoil~=2 & airfoil~=3 & airfoil~=4 & airfoil~=5 
dispC ') 

dispC *** Enter a 1, 2, 3, 4, or 5 ***') 
else 
flag=0; 
end 
end 
end 

thetao=input('collective pitch at .7 r/R (deg): ')/57.3; 
while isempty(thetao), 
dispC ') 

dispCYou must enter a numerical value') 
thetao=input('Collective pitch at .7 r/R (deg): ')/57.3; 
end 

Afh=input('Aircraft equivalent flatplate area (fC2):'); 
while isempty(Afh) 
dispC ') 

dispCYou must enter a numerical value') 

Afh=input('Aircraft equivalent flatplate area (ft^2):'); 
end 

Afv=input('Vertical projected area (fC2): '); 
while isempty(Afv) 
dispC ') 

dispCTou must enter a numerical value') 
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Afv=input('Vertical projected area (ft^2): '); 
end 

Shoriz=input('Horizontal tail area, 0 if none (ft^2):'); 
while isempty(Shoriz), 
dispC ') 

dispCYou must enter a numerical value') 
Shoriz=input('Horizontal tail area, 0 if none (ft^2): '); 
end 

if Shoriz~=0, 

bhoriz=input('Horizontal tail span (ft): '); 
while isempty(bhoriz), 
dispC ') 

dispCYou must enter a numerical value') 
bhoriz=input('Horizontal tail span (ft): '); 
end 

if bhoriz < le-10, 
bhoriz=le-10; 
end 

CLhoriz=input('Expected CL for the horizontal tail: '); 
while isempty(CLhoriz), 
dispC ') 

dispCYou must enter a numerical value') 
CLhoriz=input('Expected CL for the horizontal tail: '); 
end 

CDohoriz=input('Horizontal tail profile drag coef (CDo): '); 
while isempty(CDohoriz), 
dispC ') 

dispCYou must enter a numerical value') 
CDohoriz=input('Horizontal tail profile drag coef (CDo): '); 
end 
else 

Shoriz= le-10; 
bhoriz= le-10; 

CLhoriz=0; 

CDohoriz=0; 

end 

Svert=input('Vertical tail area, 0 if none (fC2): '); 
while isempty(Svert), 
dispC ') 

dispCYou must enter a numerical value') 

Svert=input('Vertical tail area, 0 if none (ft''2): '); 
end 

if Svert~=0 

bvert=input('Vertical tail span (ft): '); 
while isempty(bvert), 
dispC ') 

dispCYou must enter a numerical value') 
bvert=input('Vertical tail span (ft): '); 
end 

if bvert < le-10, 
bvert= le-10; 
end 

CLvert=input('Expected CL for the vertical tail: '); 
while isempty(CLvert), 
dispC ') 

dispCYou must enter a numerical value') 
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CLvert=input('Expected CL for the vertical tail: '); 
end 

CDovert=input('Vertical tail profile drag coef (CDo):'); 
while isempty(CDovert), 
dispC ') 

dispCYou must enter a numerical value') 
CDovert=input('Vertical tail profile drag coef (CDo):'); 
end 
else 

Svert=le-10; 

bvert=le-10; 

CLvert=0; 

CDovert=0; 

end 

twoPinput=input('2P longitudinal input (degs): '); 
while isempty(twoPinput), 
dispC ') 

dispCYou must enter a numerical value') 
twoPinput=input('2P longitudinal input (degs): '); 
end 

threePinput=input('3P input (degs): '); 
while isempty(threePinput), 
dispC ') 

dispCYou must enter a numerical value') 
threePinput=input('3P input (degs): '); 
end 

threePoffset=input('3P offset, must be between 0 and 119 degs: '); 
while isempty(threePoffset) or threePoffset>119 or threePoffset<0, 
dispC) 

dispCYou must enter a numerical value between 0 and 119 degs') 
threePoffset=input('3P offset (degs): '); 
end 

Swing=le-10; 

bwing=le-10; 

CLwing=0; 

CDowing=0; 

ewing=le-10; 

Taux=0; 

rvrCP=2; 

if helochoice==2 

Taux=input('Auxiliary thrust (lbs): '); 
while isempty(Taux), 
dispC ') 

dispCYou must enter a numerical value') 

Taux=input('Auxiliary thrust (lbs): '); 
end 

Swing=input('Wing area, 0 if no wing (ft^2): '); 
while isempty(Swing) 
dispC ') 

dispCYou must enter a numerical value') 

Swing=input('Wing area, 0 if no wing (fC2): '); 
end 

if Swing~=0, 

bwing=input('Wing span (ft):'); 
while isempty(bwing), 
dispC ') 
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dispCYou must enter a numerical value') 
bwing=input('Wing span (ft):'); 
end 

if bwing < le-10, 
bwing=le-10; 
end 

CLwing=input('Expected CL for the wing: '); 
while isempty(CLwing), 
dispC ') 

dispCYou must enter a numerical value') 

CLwing=input('Expected CL for the wing: '); 
end 

CDowing=input('Wing profile drag coef (CDo): '); 
while isempty(CDowing), 
dispC ') 

dispCYou must enter a numerical value') 

CDowing=inputCWing profile drag coef (CDo): '); 
end 

ewing=input('Wing efficiency factor (e): '); 
while isempty(ewing), 
dispC ') 

dispCYou must enter a numerical value') 
ewing=input('Wing efficiency factor (e): '); 
end 

if ewing < le-10, 
ewing= le-10; 
end 
else 

Swing= le-10; 
bwing= le-10; 

CLwing=0; 

CDowing=0; 

ewing=le-10; 

end 

GW1=GW; 

clc 

ifVinf> 16.9 
dispC) 

dispCDo you want to fix the Tip Path Plane (recommended for compound 

helicopters),') 

alphaTfix=input('l. yes or 2. no?'); 
if alphaTfix==l 
dispC) 

alphaTinput=input('Enter the Tip Path Plane angle in degrees: '); 
alphaTinput=alphaTinput/57.3; 

dispC*** You have fixed the Tip Path Plane angle, please set auxiliary ***') 
dispC*** thrust to equal total drag to avoid erroneous solutions. ***') 
else 
dispC) 

dispCTip path plane will be determined by JANRAD') 
dispCPress any key to continue...') 
pause 
end 

% added 14-16 Sep 2003 to allow user to disable the retrim routine, select Taux based 
% on drag, and allow for RVR based cyclic inputs in trim.m 
flag=l; 
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while flag > 0 
dispC) 

dispCDo you want JANRAD to set auxiliary thrust to equal total drag,') 
Tauxchoice=input('l. yes or 2. no?'); 
while isempty(Tauxchoice), 
dispC ') 

dispCYou must enter a 1 or 2') 

dispCDo you want JANRAD to set auxiliary thrust to equal total drag,') 
Tauxchoice=input('l. yes or 2. no?'); 
end 

if Tauxchoice~=l & Tauxchoice~=2 
dispC ') 

dispC *** Enter a 1 or 2 ***') 
elseif Tauxchoice==l 

dispCJANRAD will adjust auxiliary thrust to equal total drag') 

dispCPress any key to continue...') 

pause 

flag=0; 

elseif Tauxchoice==2 

dispC*** Auxiliary thrust defined by user. ***') 
dispCPress any key to continue...') 
pause 
flag=0; 
else 
flag=0; 
end 
end 

rvrCP=2; 

ifVinf/1.69>200 

rvrCP=0; 

flag=l; 

while flag > 0 
dispC) 

dispCDo you want JANRAD to use RVR based (constrained) 2/rev Cyclic Pitch, 

') 

rvrCP=inputCl. yes or 2. no?'); 
while isempty(rvrCP), 
dispC') 

dispCYou must enter a 1 or 2') 

dispCDo you want JANRAD to use RVR based (constrained) 2/rev Cyclic 

Pitch,') 

rvrCP=input('l. yes or 2. no?'); 
end 

if rvrCP~=l & rvrCP~=2 
dispC') 

dispC *** Enter a 1 or 2 ***') 
elseif rvrCP==l 
dispC) 

dispCRVR 2/rev Cyclic Pitch will be used') 
dispCPress any key to continue...') 
pause 
flag=0; 

elseif rvrCP==2 
dispC) 

dispCConventional 2/rev cyclic pitch will be used') 
dispCPress any key to continue...') 
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pause 

flag=0; 

else 

flag=0; 

end 

end 

else 

end 

else 

end 

else 

end 

clc 

dispC ') 
dispC ') 

dispC *** Save Instructions ***') 

dispC ') 

dispC A. Save the data to a specified file name.') 

dispC B- Do not use an extension or qoutations.') 

dispC C. Use letter/number combinations of 6 characters or less.') 

dispC D. The file will be saved with a ".mat" extension.') 

dispC ') 

dispC ex: dsgn_r) 
dispC ') 

dispC E. If you do not enter a name, the default is "dsgn_r") 

flag=l; 

while flag > 0 

filenamel=inputC save file as: ','s'); 
if isempty(filenamel) 
filename 1='dsgn_r; 
end 

if length(filenamel) > 20, 
dispC') 

dispCUse 20 characters or less.') 
flag=l; 
else 
flag=0; 
end 
end 

eval(['save ',filenamel]) 
check=0; 
else 
dispC ') 

dispC Enter a 1 or 2') 
end 
end 
clc 

dispC ') 
pause(3) 
check2=l; 
while check2 > 0 
clc 

dispC') 

dispC *** EXECUTION MENU ***') 

dispC') 

dispC 1 • Rotor Performance Analysis') 
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dispC 2. Change Data') 

dispC 3. Quit') 

check3=l; 
while check3 > 0 

answer=input(' Enter a 1, 2, or 3 »'); 
save temp check 1 check2 check3 filename 1 
if answer==l, 
perf 

load temp 
check3=0; 
elseif answer==2 
clc 

check2=0; 
check3=0; 
elseif answer==3, 
dispC ') 

dispC Thank you for using JANRAD') 
checkl=0; 
check2=0; 
check3=0; 
end 
end 
end 
end 
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APPENDIX B. DMCALC.M 


% dmcalc calculates the total drag along a blade at 
% each azimuth (psi) location 

Up=zeros(size(psi *r)); 

Ut=zeros(size(Up)); 

alpham=zeros(size(Up)); 

% added for plotting routine 
alphamplot=zeros(size(Ut)); 

dD=zeros(size(Up)); 

ddD=zeros(size(psi)); 

ddDM=zeros(size(psi)); 

for i=l;length(psi), 

Up(i,:)=vi.*cos(betao)+Vinf*sin(alphaT)*cos(betao)+Vinf*cos(alphaT)*sin(betao)*cos(psi(i)); 

Ut(i,:)=r.*omega+Vinf*cos(alphaT)*sin(psi(i)); 

mach = (Vtip.*cos(alphaT).*r./R+Vinf.*sin(psi(i)))/(49.05*sqrt(temp+460)); 
phi=atan2(Up(i,:),Ut(i,:)); 

% Replace IP theta with 2P theta 
% alpha=theta(i)+betat-phi; 
alpha=thetaXP(i)+betat-phi; 

% Test for plotting routine 
alpha 1 =thetaXP(i)+betat-phi; 
alphamplot(i,;)=alphal; 

alpham(i,:)=alpha; 


if airfoil==l, 

[CL,CD]=hh02clcd(alpha); 

% additional airfoils added for experimentation 
elseif airfoil==2, 

[CL,CD]=vrl2clcd(alpha); 
elseif airfoil==3, 

[CL,CD]=scl095r8clcd(alpha,mach); 
elseif airfoil==4, 

[CL,CD] =oo 12clcd(alpha,mach); 
elseif airfoil==5, 

[CL,CD] =rvrclcd(alpha,mach); 
end 

% Populate CLplot and CDplot matrix 
CLplot(i,:)=CL; 

CDplot(i,:)=CD; 

dD(i,:)=0.5*rho*cblade*dr*(Up(i,:).^2).*(CL.*sin(phi)+CD.*cos(phi)); 

dDM=dD(i,:).*r; 

DMpsi(i)=sum(dDM); 

% *** calculations for tip loss area *** 
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Uptip=Vinf*sin(alphaT)*cos(betao)+Vinf*cos(alphaT)*sin(betao)*cos(psi(i)); 

Uttip=(R-(r-RefQ/2)*omega+Vinf*cos(alphaT)*sin(psi(i)); 

phitip=atan2(Uptip,Uttip); 

ddD(i)=0.5*rho*cblade*(R-RefQ*(Uptip^2+Uttip.^2)*(.009*cos(phitip')); 

ddDM(i)=ddD(i)*(R-(R-Reff)/2); 

DMpsi(i)=DMpsi(i)+ddDM(i); 

end 
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APPENDIX C. PERF.M 


% *** Performance output subroutine *** 

% *** Clearing all variable in the MATLAB environment *** 
clear 

% *** Loading input and check parameters *** 
load temp 

eval(['load ',filenamel]) 
trim 

dispCtrimming complete') 

% *** Calculation of output parameters *** 

Protor=mean(DMpsi)*b*omega/550; 

Qrotor=mean(DMpsi)*b; 
solidity=b * cblade/(pi * R); 

CQ=Qrotor/(Adisk*rho*Vtip^2*R); 

CH=Hrotor/(Adisk*rho*Vtip''2); 

CT_sig=CT/ solidity; 

CQ_sig=CQ/solidity; 

CH_sig=CH/solidity; 

Machtip=( V tip*cos(alphaT)+Vinf)/(49.05 * sqrt(temp+460)); 
ifVinf< 16.9, 

DL=T/(pi*R^2); 

FM=(T*sqrt(DL/(2*rho)))/(550*Protor); 

else 

DL=0; 

FM=0; 

end 

clc 

dispC ') 

dispC * * * PERFORMANCE CALCULATIONS COMPLETE * 

dispC ') 

dispC Do you want the performance results displayed on screen?') 

flag=l; 

while flag > 0 

answer=input(' 1. yes 2. no »'); 
if answer==2, 
flag=0; 

elseif answer==l 
% *** output to screen *** 
clc 

dispC ') 

dispC *** INPUT DATA ***') 
eval(['disp(" ',filenamel,'")']) 
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dispC ') 

fprintfC Forward velocity = %6.0f kts\n',Vinf/1.69) 
fprintfC Temperature = %6.0f degs F\n',temp) 

fprintfC Pressure Altitude = %6.0f ft\n',PA) 
fprintfC Gross weight = %6.0f lbs\n',GW) 

fprintfC Number of blades = %6.0f \n',b) 
fprintfC Rotor radius = %6.2f ft\n',R) 

fprintfC Blade chord = %6.2f ft\n',cblade) 

fprintfC Blade twist = %6.2f degs\n',-l*twist*57.3) 

if airfoil==l 

dispC Blade airfoil = HH-02') 

elseif airfoil==2 

dispC Blade airfoil = VR-12') 

elseif airfoil==3 

dispC Blade airfoil = SC1095') 

elseif airfoil==4 

dispC Blade airfoil = 0012') 

elseif airfoil==5 

dispC Blade airfoil = RVR') 

end 

fprintfC Blade lift curve slope = %6.2f \n',a) 

fprintfC Blade wieght = %6.2f lbs\n',wblade) 

fprintfC Rotational velocity = %6.2f rads/sec\n',omega) 

fprintfC Blade grip length = %6.2f ft\n',grip) 

fprintfC Hinge offset = %6.2f ft\n',e) 

dispC ') 

dispC ') 

dispCpress any key to continue...') 
pause 

clc 

dispC ') 

dispC *** INPUT DATA CONTINUED ***') 
eval(['disp(" ',filenamel,"')']) 

dispC ') 

fprintfCEquivalent flat plate area = %6.2f ft''2\n',Afh) 
fprintfC Vertical Projected area = %6.2f fU2\n',Afv) 
fprintfC Wing Area = %6.2f ft^2\n',Swing) 

fprintfC Wing Span = %6.2f ft\n',bwing) 

fprintfC Wing CL = %6.2f \n',CLwing) 

fprintfC Wing CDo = %6.4f \n',CDowing) 

fprintfC Wing efficiency factor = %6.2f \n',ewing) 
fprintfC Horizontal tail area = %6.2f ft''2\n',Shoriz) 

fprintfC Horizontal tail span = %6.2f ft\n',bhoriz) 

fprintfC Horizontal tail CL = %6.2f \n',CLhoriz) 

fprintfC Horizontal tail CDo = %6.4f \n',CDohoriz) 
fprintfC Vertical tail area = %6.2f fU2\n',Svert) 
fprintfC Vertical tail span = %6.2f ft\n',bvert) 
fprintfC Vertical tail CL = %6.2f \n',CLvert) 
fprintfC Vertical tail CDo = %6.4f \n',CDovert) 

fprintfC Auxiliary thrust = %6.0f lbs\n',Taux) 
fprintfC Vertical Auxiliary thrust = %6.0f lbs\n',vertaux) 
fprintfC2/rev cyclic input (deg) = %6.0f \n',theta2c*57.3) 
fprintfC3/rev cyclic input (deg) = %6.2f \n',theta3s*57.3) 
fprintf('3/rev offset (deg) = %6.2f \n',theta3Poffset*57.3) 
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dispC ') 
dispC ') 

dispCpress any key to continue...') 
pause 

clc 

dispC ') 

dispC *** CALCULATED DATA ***') 
eval(['disp(" ',filenamel,"')']) 

dispC ') 

fprintfC Fuselage drag = %6.0f lbs\n',Dfuse) 
fprintfC Rotor drag = %6.0f lbs\n',Hrotor) 

fprintfC Wing Lift = %6.0f lbs\n',Lwing) 

fprintfC Wing Drag = %6.0f lbs\n',Dwing) 

fprintfC Horizontal tail lift = %6.0f lbs\n',Lhoriz) 
fprintfC Horizontal tail drag = %6.0f lbs\n',Dhoriz) 
fprintfC Vertical tail side force = %6.0f lbs\n',Lvert) 
fprintfC Vertical tail drag = %6.0f lbs\n',Dvert) 
fprintfC Tip path angle = %6.2f degs\n',alphaT*57.3) 
fprintfC Rotor Coning angle = %6.2f degs\n',betao*57.3) 
fprintfCLocation of mean thrust(r/R)= %6.2f \n',rT2) 
fprintfCCollective pitch at .7 r/R = %6.2f degs\n',thetao*57.3) 
fprintfClst lat cyclic term-Al(deg)= %6.2f \n',thetalc*57.3) 
fprintfClst long cyclic term-B 1 (deg)= %6.2f \n',thetals*57.3) 
dispC ') 
dispC ') 

dispCpress any key to continue...') 
pause 

clc 

dispC ') 

dispC *** CALCULATED DATA CONTINUED ***') 
eval(['disp(" ',filenamel,'")']) 

dispC ') 

fprintfC solidity (sigma) = %6.3f\n',solidity) 

fprintfC Disk Loading = %6.2f lbs/ft^2\n',DL) 

fprintfC Figure of Merit = %6.2f \n',FM) 
fprintfC CT/sigma = %6.3f \n',CT_sig) 

fprintfC CQ/sigma = %6.4f \n',CQ_sig) 

fprintfC CH/sigma = %6.4f \n',CH_sig) 

fprintfCTip Mach number of adv blade= %6.3f \n',Machtip) 
fprintfC Advance ratio = %6.3f \n',mu) 

fprintfCRotor thrust required (TPP)= %6.0f lbs\n',T) 
fprintfC Rotor power required= %6.0f h.p.\n',Protor) 
fprintfC Rotor torque= %6.0f ft-lbs\n',Qrotor) 

dispC ') 
dispC ') 

dispCpress any key to continue...') 
pause 

clc 

flag=0; 
else 
dispC ') 

dispC Enter a 1 or 2') 
end 
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end 


% *** output to disk (text file) *** 
dispC ') 

dispC *** Saving data to disk ***') 

pause(2) 

diary on 

diary off 

delete diary 

diary on 

dispC *** RESULTS ***') 

eval(['disp(" ',filenamel,"')']) 

fprintfC Forward Velocity = %6.0f kts\n',Vinf/1.69) 

fprintfC Temperature = %6.0f degs F\n',temp) 

fprintfC Pressure Altitude = %6.0f ft\n',PA) 
fprintfC Gross weight = %6.0f lbs\n',GW) 

fprintfC Munber of blades = %6.0f \n',b) 

fprintfC Rotor Radius = %6.2f ft\n',R) 

fprintfC Blade Chord = %6.2f ft\n',cblade) 

fprintfC Blade twist = %6.2f degs\n',-l*twist*57.3) 

fprintfC Blade lift curve slope = %6.2f \n',a) 
fprintfC Blade weight = %6.2f lbs\n',wblade) 

fprintfC Rotational Velocity = %6.2f rads/sec\n',omega) 
fprintfC Blade grip length = %6.2f ft\n',grip) 
fprintfC Hinge offset = %6.2f ft\n',e) 
fprintfC Equivalent flat plate area= %6.2f ft''2\n',Afh) 
fprintfC Vertical projected area = %6.2f fU2\n',Afv) 
fprintfC Wing Area = %6.2f ft''2\n',Swing) 

fprintfC Wing Span = %6.2f ft\n',bwing) 

fprintfC Wing CL = %6.2f \n',CL wing) 

fprintfC Wing CDo = %6.4f \n',CDowing) 

fprintfC Wing Efficicency Factor = %6.2f \n',ewing) 
fprintfC Horizontal tail area = %6.2f ft^2\n',Shoriz) 
fprintfC Horizontal tail span = %6.2f ft\n',bhoriz) 
fprintfC Horizontal tail CL = %6.2f \n',CLhoriz) 
fprintfC Horizontal tail CDo = %6.4f \n',CDohoriz) 
fprintfC Vertical tail area = %6.2f fU2\n',Svert) 

fprintfC Vertical tail span = %6.2f ft\n',bvert) 

fprintfC Vertical tail CL = %6.2f \n',CLvert) 

fprintfC Vertical tail CDo = %6.4f \n',CDovert) 
fprintfC Fuselage drag = %6.0f lbs\n',Dfuse) 
fprintfC Rotor drag = %6.0f lbs\n',Drotor) 

fprintfC Wing lift = %6.0f lbs\n',Lwing) 

fprintfC Wing drag = %6.0f lbs\n',Dwing) 

fprintfC Horizontal tail lift = %6.0f lbs\n',Lhoriz) 

fprintfC Horizontal tail drag = %6.0f lbs\n',Dhoriz) 
fprintfC Vertical tail side force= %6.0f lbs\n',Lvert) 
fprintfC Vertical tail drag = %6.0f lbs\n',Dvert) 
fprintfC Auxiliary thrust = %6.0f lbs\n',Taux) 
fprintfC Vertical Auxiliary thrust = %6.0f lbs\n',vertaux) 
fprintfC Tip path angle = %6.2f degs\n',alphaT*57.3) 
fprintfC Rotor coning angle = %6.2f degs\n',betao*57.3) 
fprintfCLocation of mean thrust(r/R)= %6.2f \n',rT2) 
fprintfC Collective pitch at .7 r/R= %6.2f degs\n',thetao*57.3) 
fprintf(Tst lat cyclic term-Al(deg)= %6.2f \n',thetalc*57.3) 
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fprintf('lst long cyclic term-B 1 (deg)= %6.2f \n',thetals*57.3) 
fprintf('2/rev cyclic term (deg) = %6.2f \n',theta2c*57.3) 
fprintf('3/rev cyclic term (deg) = %6.2f \n',theta3s*57.3) 
fprintf('3/rev offset (deg) = %6.2f \n',theta3Poffset*57.3) 
fprintfC solidity = %6.3f\n',solidity) 

fprintfC Disk loading = %6.2f lbs/ft^2\n',DL) 

fprintfC Figure of Merit = %6.2f \n',FM) 
fprintfC CT/sigma = %6.3f \n',CT_sig) 

fprintfC CQ/sigma = %6.4f \n',CQ_sig) 

fprintfC CH/sigma = %6.4f \n',CH_sig) 

fprintfCTip Mach of the adv. blade = %6.3f \n',Machtip) 
fprintfC Advance Ratio = %6.3f \n',mu) 
fprintfCRotor thrust required(TPP) = %6.0f lbs\n',T) 
fprintfC Rotor power required = %6.0f h.p.\n',Protor) 
fprintfC Rotor Torque = %6.0f ft-lbs\n',Qrotor) 
diary off 

% *** Output to disk (.mat file containing matrix variables: 

% r, psi, vi, theta, thetaXP, betat, alpha, Tpsi, Mpsi, DMpsi, dT, dM 

% dD, CLplot, CDplot) *** 

% *** Configuring variables for output *** 

theta=theta*57.3; 

thetaXP=thetaXP*57.3; 

betat=[betattwist*(0.7-(Reff+(R-Reff)/2)/R)]*57.3; 

alpha=alpham*57.3;,alpha=[alpha zeros(size(psi))]; 

Mpsi=Mpsi(:,length(Mpsi(l,:))-!); 

dM=[dM ddM]; 

psi=psi*57.3; 

r=[r (R-(R-Reff)/2)]; 

vi=[vi 0]; 

dr=dr*12; 


clc 

flagl=l; 
while flagl > 0 
dispC ') 

wantplots=inputCDo you want to view plots, 1. yes or 2. no?'); 
if wantplots==2, 
flag 1=0; 

elseif wantplots==l, 
plotter 
flag 1=0; 
else 
dispC ') 

dispC Enter a 1 or 2') 
end 
end 

clc 

dispC ') 

dispC * * * PERFORMANCE OUTPUT DATA INSTRUCTIONS * * * 

dispC ') 

dispC A. Single value data saved to a file named:') 
eval(['disp(" '",filenamel,'.prf"")']) 


93 



dispC ') 

dispC B. This is a text file, use the tyoe command to view the') 
dispC file or use a text editor to view/print the file.') 
eval( ['flag=exist('",filename 1,' .prf);']) 
if flag < 1, 

eval(['!rename diary ',filenamel,'.prf"]) 
else 

eval(['!del ',filenamel,'.prf"]) 
eval(['!rename diary ',filenamel,'.prf"]) 
end 

filename2= [filename 1 '_p']; 
dispC ') 

dispC C. Matrix and vector data saved to a file named:') 
eval(['disp(" '”,filename2,'.mat"")']) 
dispC ') 

dispC D. This is a ".mat" binary file, use the load command') 
dispC to retrieve the data for plotting.') 

eval(['save ',filename2,' r psi vi theta thetaXP betat alpha Tpsi Mpsi DMpsi dT dM dD CLplot 

CDplof]) 

dispC ') 

dispC *** end of performance ***') 
dispC ') 
dispC ') 

dispC press any key to continue...') 
pause 
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APPENDIX D. PLOTTER.M 


% *** Plotting Routine added 28 July 2003 some code copied/modified from lines 28-72 of 

create_plots 

% JANRAD v6.0 *** 
clc 

dispC ') 

dispC * * * * Available Plots ****') 

dispC ') 

dispC 1 • Rotor Azimuth Angle v. Cyclic Pitch') 
dispC 2. Blade Azimuth Position v. Airload') 
dispC 3. Rotor Airload Distribution') 
dispC 4. Blade Azimuth Position v. Lift') 
dispC 5. Rotor Lift Distribution') 

dispC 6. Angle of Attack Distribution') 

dispC 7. Rotor Tip Path Plane Stall pattern and Thrust Distribution') 

dispC 8. Rotor Trim Lissajous Figure') 

dispC 9. Lift Distribution (Contour Plot)') 

dispC 10. CL Distribution (Contour Plot)') 

dispC 11 • Mach Number Distribution (Contour Plot)') 

dispC ') 

dispC ') 

flag2=l; 

drinches=dr; 

while flag2 > 0 

plotchoice=input('Please choose a plot by entering 1-11.'); 

if plotchoice== 1, 
flag2=0; 

if twoPinput==0 & threePinput==0 
figure(l) 
plot(psi,thetaXP) 
axis([0,360,-10,30]); 
legend)'IP Harmonic Cyclic Pitch') 

title(['Rotor Azimuth Angle v. Harmonic Cyclic Pitch for ',num2str(Vinf/L68781),... 

' Kts, advance ratio=',num2str(mu)]) 
xlabel('Rotor Azimuth Angle (degs)'),ylabel('Cyclic Pitch (degs)') 
elseif rvrCP==l 
figure) 1) 

thetao=thetao.*ones(size(theta2Pr)); 

plot(psi,thetao*57.3,'-',psi,(-theta-thetao*57.3),':',psi,theta2Pr*57.3,'—',psi,thetaXP,'*') 
axis([0,360,-25,20]); 

legend('Collective Pitch','IP Cyclic Pitch','2P Cyclic Pitch','Coll + IP H- 2P signal') 
title(['Rotor Azimuth Angle v. RVR Harmonic Cyclic Pitch for ',num2str(Vinf/L68781),... 

' Kts, advance ratio=',num2str(mu),', 2P input=',num2str(twoPinput),... 

' degs']) 

xlabel('Rotor Azimuth Angle (degs)'),ylabel('Cyclic Pitch (degs)') 
else 

figure) 1) 

plot(psi,theta,psi,theta2Pr*57.3,'—',psi,theta3Pr*57.3,':',psi,thetaXP,'*') 
axis([0,360,-15,35]); 

legend('lP Cyclic Pitch','2P Cyclic Pitch',... 

'3P Cyclic Pitch','lP H- 2P H- 3P signal') 

title(['Rotor Azimuth Angle v. Harmonic Cyclic Pitch for ',num2str(Vinf/L68781),... 
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' Kts, advance ratio=',num2str(mu),', 2P input=',num2str(twoPinput),... 

' degs, 3P input=',num2str(threePinput),... 

' degs, 3P offset=',num2str(threePoffset),' degs']) 
xlabel('Rotor Azimuth Angle (degs)'),ylabel('Cyclic Pitch (degs)') 
end 

elseif plotchoice==2, 
flag2=0; 
figure(2) 

subplot(3,3,8) 

plot(r./R,dT(l,:)./drinches,'k'),grid 
title ('Psi = 0 deg ') 

xlabel('Blade Position r/R');ylabel('Airload (Lb/in)'); 
subplot(3,3,2) 

plot(r./R,dT(floor(naz/2),:)./drinches,'k'),grid 
title('Psi =180 deg') 

xlabel('Blade Position r/R'); ylabel('Airload (Lb/in)'); 
subplot(3,3,6) 

plot(r./R,dT(floor(naz/4),;)./drinches,'k'),grid 
title('Psi = 90 deg ') 

xlabel('Blade Position r/R');ylabel('Airload (Lb/in)'); 
subplot(3,3,4) 

plot(r./R,dT(floor(3*naz/4),:)./drinches,'k'),grid 
title('Psi = 270 deg ') 

xlabel('Blade Position r/R'); ylabel('Airload (Lb/in)'); 

gtext(['Blade Position v. Airload for ',num2str(Vinf/1.68781),' Kts, advance 
ratio=',num2str(mu),... 

', 2P input=',num2str(twoPinput),' degs, 3P input=',... 
num2str(threePinput),' degs, 3P offset=',num2str(threePoffset),' degs']) 

elseif plotchoice==3, 
flag2=0; 
figure(3) 

[th,rl] = meshgrid((-180:360/naz:180)*pi/180,r/R); 

[x,y] = pol2cart(th,rl); 
dTl=[dT; dT(l,:)]; 
for i=l:naz+l 

dTl(i,:)=dTl(i,:)./(drinches); 

end 

mesh(x',y',dTl) 
view(315,60) 

xlabel('Starboard'); ylabel('Aft'); zlabel('Aero Load, Lb/in'); 

title(['Airload Distribution At ',num2str(Vinf/l .68781),' Kts, advance ratio=',num2str(mu),... 

', 2P input=',num2str(twoPinput),' degs, 3P input=',... 
num2str(threePinput),' degs, 3P offset=',num2str(threePoffset),' degs']) 

elseif plotchoice==4, 
flag2=0; 
figure(4) 
subplot(3,3,8) 
plot(r./R,dN(l,:)),grid 
title('Psi = 0 deg ') 
xlabel('r/R');ylabel('Lift,Lb'); 
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subplot(3,3,2) 

plot(r./R,dN(floor(naz/2),:)),grid 
title('Psi =180 deg ') 
xlabel('r/R'); ylabel('Lift,Lb'); 

subplot(3,3,6) 

plot(r./R,dN(floor(naz/4),:)),grid 
title('Psi = 90 deg ') 
xlabel('r/R');ylabel('Lift,Lb'); 

subplot(3,3,4) 

plot(r./R,dN(floor(3 *naz/4), :)),grid 
title('Psi = 270 deg ') 
xlabel('r/R'); ylabel('Lift,lb'); 

gtext(['Blade Position v. Lift for',num2str(Vinf/l.68781),' Kts, advance ratio=',num2str(mu),... 
', 2P input=',num2str(twoPinput),' degs, 3P input=',... 
num2str(threePinput),' degs, 3P offset=',num2str(threePoffset),' degs']) 

elseif plotchoice==5, 
flag2=0; 
figure(5) 

[th,rl] = meshgrid((-180:360/naz:180)*pi/180,r/R); 

[x,y] = pol2cart(th,rl); 
dNl=[dN; dN(l,:)]; 
for i=l:naz+l 
dNl(i,:)=dNl(i,:); 
end 

mesh(x',y',dNl) 
view(315,60) 

xlabel('Starboard'); ylabel('Aft'); zlabel('Lift, Lb'); 

title(['Lift Distribution At ',num2str(Vinf/L68781),' Kts, advance ratio=',num2str(mu),... 

', 2P input=',num2str(twoPinput),' degs, 3P input=',... 
num2str(threePinput),' degs, 3P offset=',num2str(threePoffset),' degs']) 

elseif plotchoice==6, 
flag2=0; 
figure(6) 

alpha=alpham*57.3; 

alpha=alpha'; 

drl=(R-e)/nbe; 

rl=e:drl;R-drl; 

rl=rl+drl/2; 

[th,r2]=meshgrid((0:360/naz:360)*pi/180,rl/R); 

[x,y]=pol2cart(th,r2); 

alpha=[alpha,alpha(:,l)]; 

% The vector v may be altered to show user defined contour lines 
ifrvrCP==l 

v=[-196,-194,-190,-180,-170,-18,-14,-10,-8,-4,-2,0,2,4,8,10,14,18]; 

else 

v=[18,-16,-14,-12,-10,-8,-6,-4,-2,0,2,4,6,8,10,12,14,16]; 

end 

hl=polarl([0,2*pi], [0,1]); 
delete (hi) 
hold on 
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hub area 


% Either [C,h2] equation may be used. The second equation, shown in line 162, depicts the 


% [C,h2]=contour(x,y,alpha,v); 

[C,h2]=contour(x(round(grip* 12/dr):nbe, 1 ;nbe+ l),y(round(grip* 12/dr):nbe, 1 :nbe+ l),alpha(round(grip* 12/ 
dr):nbe,l;nbe+l),v); 

clabel(C,h2); 

title(['Angle of Attack Distribution At ',num2str(Vinf/l.68781),' Kts, advance 
ratio=',num2str(mu),... 

', 2P input=',num2str(twoPinput),' degs, TPP angle=',... 

num2str(alphaT*57.3),' degs, 3P input=',num2str(threePinput),' degs, 3P offset=',... 
num2str(threePoffset),' degs']) 
hold on 

thetahub=[0:pi/270:2*pi]; 
hub=(grip/R)./(l+0*cos(thetahub)); 
for i=l:.1:20 

polar(thetahub,hub/i); 
hold on 
end 

hold off 
view(90,90) 

elseif plotchoice==7, 
flag2=0; 
figure(7) 

subplot(2,l,l) 

plot(psi,alpha(naz,l:naz)),set(gca,'XTick', [0:90:360],'YTick', [-8:2:14]) 
xlabel('Azimuth Angle, degs'),ylabel('Tip Angle of Attack, degs') 
title(['Rotor Tip Path Plane Stall Pattern, advance ratio=',num2str(mu)]) 
text(120,10,'Ideal') 
hold on 

psiplotl= [0,45,45,135,135,180,210,250,270,290,360]; 
ideal= [10,7,-.5,-.5,7,10,ll,12,12.5,13,10.5]; 
plot(psiplotl,ideal,'r—') 
hold off 


dTplot=sum(dT')/(rho*Adisk*(omega*R)''2); 
subplot(2,1,2) 

plot(psi,dTplot),set(gca,'XTick', [0:90:360]) 

xlabel('Azimuth Angle, degs'),ylabel('Thrust Coefficient (one blade)') 
title(['Actual and Ideal Thrust Distribution at Stall, advance ratio=',num2str(mu)]) 
text(180,.003,'Ideal') 
hold on 

psiplot2= [0 45 45 135 135 180 210 240 270 315 340 360]; 

ideal2= [.003 .0048 .0001 .0001 .0048 .003 .002 .001 .0009 .0015 .002 .0029]; 
plot(psiplot2,ideal2,'r—') 
hold off 

gtext([num2str(Vinf/l.68781),' Kts, advance ratio=',num2str(mu),... 

', 2P input=',num2str(twoPinput),' degs, TPP angle=',... 

num2str(alphaT),' degs, 3P input=',num2str(threePinput),' degs, 3P offset=',... 

num2str(threePoffset),' degs']) 

elseif plotchoice==8, 
flag2=0; 
figure(8) 
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Mpsil=[Mpsi; Mpsi(l,l)]; 

% if desired convert to in * lbs 

Mpsil=Mpsil*12; 

theta 1=[theta; theta(l,l)]; 

plot(Mpsi 1, theta 1) 

axis([-le6,3e6,0,35]); 

theta 1 =theta 1 *pi/l 80; 

alphaT=alphaT*57.3 

% determine Lissajous area 
area=trapz(Mpsi* 12,theta/57.3); 
area=abs(area); 

xlabel('Thrust Moment, in lb'),ylabel('Blade Angle, degs') 

title(['Lissajou Figure ',num2str(Vinf/l.68781),' Kts, advance ratio=',num2str(mu),... 

', 2P input=',num2str(twoPinput),' degs, Lissajou Area=',num2str(area),' in lbs, TPP 
angle=',num2str(alphaT),' degs, 3P input=',... 

num2str(threePinput),' degs, 3P offset=',num2str(threePoffset),' degs']) 

elseif plotchoice==9, 
flag2=0; 
figure(9) 
drl=(R-e)/nbe; 
rl=e:drl;R-drl; 
rl=rl+drl/2; 

[th,r2]=meshgrid((0:360/naz:360)*pi/180,rl/R); 

[x,y] =pol2cart(th,r2); 


for i=l:naz; 
for j=l;naz; 

dT2(i,j)=dT(i,j)./(drinches); 

end 

end 

dT3=[dT2; dT2(l,:)]; 
dT3=dT3'; 

v=[-5,-4,-3,-2,-l,0,4,8,12,16,20,24,28,32,36,40,44,48,52,56,60,64,68]; 
hl=polarl([0,2*pi], [0,1]); 
delete(hl) 
hold on 

% Either [C,h2] equation may be used. The second depicts the hub area 
% [C,h2]=contour(x,y,dT3,v); 

[C,h2]=contour(x(round(grip*12/dr):nbe,l:nbe+l),y(round(grip*12/dr):nbe,l;nbe+l),dT3(round(grip*12/dr 

):nbe,l;nbe+l),v); 

clabel(C,h2); 

title(['Airload Distribution At ',num2str(Vinf/l .68781),' Kts, advance ratio=',num2str(mu),... 

', 2P input=',num2str(twoPinput),' degs, TPP angle=',... 

num2str(alphaT*57.3),' degs, 3P input=',num2str(threePinput),' degs, 3P offset=',... 
num2str(threePoffset),' degs']) 
hold on 

thetahub=[0:pi/270:2*pi]; 
hub=(grip/R)./(l+0*cos(thetahub)); 
for i=l;.1:20 

polar(thetahub,hub/i); 
hold on 
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end 

hold off 
view(90,90) 

elseif plotchoice==10, 
flag2=0; 
figure(lO) 
drl=(R-e)/nbe; 
rl=e:drl;R-drl; 
rl=rl+drl/2; 

[th,r2]=meshgrid((0:360/naz:360)*pi/180,rl/R); 

[x2,y2] =pol2cart(th,r2); 

CLplot 1 =CLplot'; 

CLplotl=[CLplotl,CLplotl(:,l)]; 

% The vector v may be altered to show user defined contour lines 

% v=[-l,-.8,-.6,-.4,-.3,-.2,-.1,0,-1,-2,.3,.4,.6,.8,1.0,1-2,1.4]; 

v=[-l,-.8,-.6,-.4,-.2,0,.2,.4,.6,.8,1.0,1.2,1.4]; 

hl=polarl([0,2*pi], [0,1]); 

delete(hl) 

hold on 

[C,h2]=contour(x2(round(grip*12/dr):nbe,l;nbeH-l),y2(round(grip*12/dr):nbe,l:nbeH-l),CLplotl(round(gri 

p*12/dr);nbe,l:nbeH-l),v); 

clabel(C,h2); 

title(['CL Distribution At ',num2str(Vinf/l.68781),' Kts, advance ratio=',num2str(mu),... 

', 2P input=',num2str(twoPinput),' degs, TPP angle=',... 

num2str(alphaT*57.3),' degs, 3P input=',num2str(threePinput),' degs, 3P offset=',... 
num2str(threePoffset),' degs']) 
hold on 

thetahub=[0:pi/270:2*pi]; 
hub=(grip/R)./(lH-0*cos(thetahub)); 
for i=l:.1:20 

polar(thetahub,hub/i); 
hold on 
end 

hold off 
view(90,90) 

elseif plotchoice== 11, 
flag2=0; 
figure(ll) 
drl=(R-e)/nbe; 
rl=e:drl:R-drl; 
rl=rlH-drl/2; 

[th,r2]=meshgrid((0:360/naz:360)*pi/180,rl/R); 

[x3 ,y 3 ] =pol2cart(th,r2); 

machplot 1 =machplot'; 

machplot2= [machplot 1 ,machplot!(:,!)]; 

v=[-.4,-.3,-.2,-.l,0,.l,.2,.3,.4,.5,.6,.7,.8,.9,l]; 
hl=polarl([0,2*pi], [0,1]); 
delete(hl) 
hold on 

[C,h2]=contour(x3,y3,machplot2,v); 
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clabel(C,h2); 

title(['Mach Number Distribution At ',num2str(Vinf/1.68781),' Kts, 
ratio=',num2str(mu)]) 
hold on 

thetahub=[0:pi/270:2*pi]; 
hub=(grip/R)./(l+0*cos(thetahub)); 
for i=l:.1:20 

polar(thetahub,hub/i); 
hold on 
end 

hold off 
view(90,90) 

else 

dispC ') 

dispC *** Enter a 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, or 11 *** ') 

end 

end 

clc 

flag3=l; 
while flag3 > 0 
dispC ’) 

wantmoreplots=input('Do you want to view additional plots, 1. yes or 2. no?'); 
if wantmoreplots==2, 
flag3=0; 

elseif wantmoreplots== 1 
plotter 
flag3=0; 
else 
dispC') 

dispC Enter a 1 or 2') 
end 
end 

% *** End of Plotting Routine *** 


advance 
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APPENDIX E. POLARl.M 


function hpol = polar(theta,rho,line_style) 

%POLAR Polar coordinate plot. 

% POLAR(THETA, RHO) makes a plot using polar coordinates of 
% the angle THETA, in radians, versus the radius RHO. 

% POLAR(THETA,RHO,S) uses the linestyle specified in string S. 
% See PLOT for a description of legal linestyles. 

% 

% See also PLOT, LOGLOG, SEMILOGX, SEMILOGY. 

% Copyright 1984-2002 The MathWorks, Inc. 

% SRevision; 5.22 $ $Date: 2002/04/08 21:44:28 $ 

if nargin < 1 

error('Requires 2 or 3 input arguments.') 
elseif nargin == 2 
if isstr(rho) 

line_style = rho; 
rho = theta; 

[mr,nr] = size(rho); 
ifmr== 1 
theta = 1 :nr; 
else 

th = (l:mr)'; 
theta = th(:,ones(l,nr)); 
end 
else 

line_style = 'auto'; 
end 

elseif nargin == 1 
line_style = 'auto'; 
rho = theta; 

[mr,nr] = size(rho); 
ifmr == 1 
theta = l:nr; 
else 

th = (l:mr)'; 
theta = th(;,ones(l,nr)); 
end 
end 

if isstr(theta) I isstr(rho) 

error('Input arguments must be numeric.'); 
end 

if ~isequal(size(theta),size(rho)) 

error('THETA and RHO must be the same size.'); 
end 

% get hold state 
cax = newplot; 

next = lower(get(cax,'NextPlot')); 
hold_state = ishold; 

% get x-axis text color so grid is in same color 

103 



tc = get(cax,'xcolor'); 

Is = get(cax,'gridlinestyle'); 


% Hold on to current Text defaults, reset them to the 

% Axes' font attributes so tick marks use them. 

fAngle = get(cax, 'DefaultTextFontAngle'); 

fName = get(cax, 'DefaultTextFontName'); 

fSize = get(cax, 'DefaultTextFontSize'); 

fWeight = get(cax, 'DefaultTextFontWeight'); 

fUnits = get(cax, 'DefaultTextUnits'); 

set(cax, 'DefaultTextFontAngle', get(cax, 'FontAngle'), ... 

'DefaultTextFontName', get(cax, 'FontName'), ... 

'DefaultTextFontSize', get(cax, 'FontSize'),... 

'DefaultTextFontWeight', get(cax, 'FontWeight'), ... 

'DefaultT extUnits', 'data') 

% only do grids if hold is off 
if ~hold_state 

% make a radial grid 
hold on; 

maxrho = max(abs(rho(:))); 

hhh=plot([-maxrho -maxrho maxrho maxrho],[-maxrho maxrho maxrho -maxrho]); 
set(gca, 'dataaspectratio', [1 1 1 ], 'plotboxaspectratiomode', 'auto') 

V = [get(cax,'xlim') get(cax,'ylim')]; 
ticks = sum(get(cax,'ytick')>=0); 
delete (hhh); 

% check radial limits and ticks 

rmin = 0; rmax = v(4); rticks = max(ticks-l,2); 
if rticks >5 % see if we can reduce the number 
if rem(rticks,2) == 0 
rticks = rticks/2; 
elseif rem(rticks,3) == 0 
rticks = rticks/3; 
end 
end 

% define a circle 
th = 0:pi/50;2*pi; 
xunit = cos(th); 
yunit = sin(th); 

% now really force points on x/y axes to lie on them exactly 
inds = l;(length(th)-l)/4;length(th); 
xunit(inds(2;2:4)) = zeros(2,l); 
yunit(inds(l;2;5)) = zeros(3,l); 

% plot background if necessary 
if ~isstr(get(cax,'color')), 
patch('xdata',xunit*rmax,'ydata',yunit*rmax, ... 
'edgecolor',tc,'facecolor',get(gca,'color'),... 

'handlevisibility','off); 

end 

% draw radial circles 

c82 = cos(82*pi/180); 
s82 = sin(82*pi/180); 
rticks= 1; 
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rinc = (rmax-rmin)/rticks; 

for i=(rmin+rinc):rinc:rmax 
hhh = plot(xunit*i,yunit*i,Is,'color',tc,'linewidth',l,... 

'handlevisibility', ’off); 
text((i+rinc/20) *c82, (i+rinc/20) * s82, ... 

[' ' num2str(i)],'verticalalignment','bottom',... 
'handlevisibility','off) 
end 

set(hhh,'linestyle','-') % Make outer circle solid 

% plot spokes 

th = (l:6)*2*pi/12; 

cst = cos(th); snt = sin(th); 

cs = [-cst; cst]; 

sn = [-snt; snt]; 

plot(rmax*cs,rmax*sn,ls,'color',tc,'linewidth',l,... 
'handlevisibility',’off) 

% annotate spokes in degrees 
rt = l.l*rmax; 
for i = l;length(th) 

text(rt*cst(i),rt*snt(i),int2str(i*30),... 

'horizontalalignment', 'center',... 

'handlevisibility','off); 
if i == length(th) 
loc = int2str(0); 
else 

loc = int2str(180H-i*30); 
end 

text(-rt*cst(i),-rt*snt(i),loc,'horizontalalignment','center',... 
'handlevisibility','off) 

end 

% set view to 2-D 
view(2); 

% set axis limits 

axis(rmax*[-l 1 -1.15 1.15]); 
end 

% Reset defaults. 

set(cax, 'DefaultTextFontAngle', fAngle , ... 
'DefaultTextFontName', fName, ... 

'DefaultTextFontSize', fSize, ... 

'DefaultTextFontWeighf, fWeight, ... 
'DefaultTextUnits',fUnits); 

% transform data to Cartesian coordinates. 

XX = rho.*cos(theta); 
yy = rho.*sin(theta); 

% plot data on top of grid 
if strcmp(line_style,'auto') 
q = plot(xx,yy); 
else 

q = plot(xx,yy,line_style); 
end 
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if nargout > 0 
hpol = q; 
end 

if ~hold_state 

set(gca,'dataaspectratio',[l 1 1]), axis off; set(cax,'NextPlot',next); 
end 

set(get(gca,'xlabel'),'visible Von') 
set(get(gca,'ylabel'),'visible','on') 
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APPENDIX F. THRCALC.M 


% thrcalc calculates the total thrast along a blade at 
% each azimuth (psi) location 
% dispCbegin thrcalc') 

Up=zeros(size(psi *r)); 

Ut=zeros(size(Up)); 
dT=zeros(size(Up)); 
ddT=zeros(size(psi)); 
dN=zeros(size(Up)); 
ddN=zeros(size(psi)); 

for i=l;length(psi), 

Up(i,:)=vi.*cos(betao)+Vinf*sin(alphaT)*cos(betao)+Vinf*cos(alphaT)*sin(betao)*cos(psi(i)); 

Ut(i,:)=r.*omega+Vinf*cos(alphaT)*sin(psi(i)); 

phi=atan2(Up(i,:),Ut(i,:)); 

mach = (Vtip.*cos(alphaT).*r./R+Vinf.*sin(psi(i)))/(49.05*sqrt(temp+460)); 

% alpha=theta(i)+betat-phi; 
alpha=thetaXP(i)+betat-phi; 


if airfoil==l, 

[CL,CD]=hh02clcd(alpha); 

% additional airfoils added for experimentation 
elseif airfoil==2, 

[CL,CD]=vrl2clcd(alpha); 
elseif airfoil==3, 

[CL,CD]=scl095r8clcd(alpha,mach); 
elseif airfoil==4, 

[CL,CD] =oo 12clcd(alpha,mach); 
elseif airfoil==5, 

[CL,CD] =rvrclcd(alpha,mach); 
end 

dT(i,:)=0.5*rho.*cblade.*dr.*(Up(i,:).''2+Ut(i,:).^2).*(CL.*cos(phi)-CD.*sin(phi)); 

Tpsi(i)=sum(dT(i,:)); 

dN(i,:)=0.5*rho.*cblade.*dr.*(Up(i,:).''2+Ut(i,:).^2).*(CL.*cos(alpha)+CD.*sin(alpha)); 

Npsi(i)=sum(dT(i,;)); 

% Populate CLplot and CDplot matrix 
CLplot(i,:)=CL; 

CDplot(i,:)=CD; 

% Populate Mach matrix 
machplot(i,;)=mach; 

% Populate alpha matrix 
alpham(i,:)=alpha; 

% *** calculations for tip loss area *** 

Uptip=Vinf*sin(alphaT)*cos(betao)+Vinf*cos(alphaT)*sin(betao)*cos(psi(i)); 

Uttip=(R-(R-Reff)/2)*omega+Vinf*cos(alphaT)*sin(psi(i)); 

phitip=atan2(Uptip,Uttip); 

alphatip=theta(i)+betat-phitip; 

%ddT(i)=0.5*rho*cblade*(R-Reff)*(Uptip^2)*(-.009*sin(phitip)); 
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ddT(i)=0.5*rho*cblade*(0.5+0.5*cos(2*psi(i)))*(R-Reff)*(Uptip^2+Uttip^2)*(- 

.009*sin(phitip)); 

Tpsi(i)=Tpsi(i)+ddT(i); 

ddN(i)=0.5*rho*cblade*(0.5+0.5*cos(2*psi(i)))*(R-Reff)*(Uptip^2+Uttip''2)*(- 

.009*sin(alphatip(i))); 

Npsi(i)=Npsi(i)+ddN(i); 

end 

machplot2=machplot'; 

CLplot2=CLplot'; 
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APPENDIX G. TMCALC.M 


% tmcalc calculates the total thrust moment along a blade 
% at each azimuth (psi) location 
% dispCbegin tmcalc') 

Up=zeros(size(psi *r)); 

Ut=zeros(size(Up)); 

dM=zeros(size(Up)); 

ddM=zeros(size(psi)); 

for i=l;length(psi), 

Up(i,;)=vi.*cos(betao)+Vinf*sin(alphaT)*cos(betao)+Vinf*cos(alphaT)*sin(betao)*cos(psi(i)); 

Ut(i,:)=r.*omega+Vinf*cos(alphaT)*sin(psi(i)); 

phi=atan2(Up(i,:),Ut(i,:)); 

% alpha=theta(i,k)+betat-phi; 
alpha=thetaXP(i,k)+betat-phi; 
if airfoil==l, 

[CL,CD]=hh02clcd(alpha); 

% additional airfoils added for experimentation 
elseif airfoil==2, 

[CL,CD]=vrl2clcd(alpha); 
elseif airfoil==3, 

[CL,CD]=scl095r8clcd(alpha,mach); 
elseif airfoil==4, 

[CL,CD] =oo 12clcd(alpha,mach); 
elseif airfoil==5, 

[CL,CD] =rvrclcd(alpha,mach); 
end 

% Populate CLplot and CDplot matrix 
CLplot(i,:)=CL; 

CDplot(i,:)=CD; 

% Populate alpham matrix 
alpham(i,:)=alpha; 


dM(i,:)=0.5*rho*cblade.*r*dr.*(Up(i,:).''2+Ut(i,:).^2).*(CL.*cos(phi)-CD.*sin(phi)); 

Mpsi(i,k)=sum(dM(i,:)); 

% *** calculations for tip loss areas *** 

Uptip=Vinf*sin(alphaT)*cos(betao)+Vinf*cos(alphaT)*sin(betao)*cos(psi(i)); 

Uttip=(R-(R-Reff)/2)*omega+Vinf*cos(alphaT)*sin(psi(i)); 

phitip=atan2(Uptip,Uttip); 

ddM(i)=0.5*rho*cblade*(R-(R-Reff)/2)*(R-Reff)*(Uptip^2+Uttip''2)*(-.009*sin(phitip)); 

Mpsi(i,k)=Mpsi(i,k)+ddM(i); 

end 

CLplot2=CLplot'; 
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APPENDIX H. TRIM.M 


% Rotor trim subroutine 


dispC ') 

dispC *** EXECUTING ROTOR TRIM ROUTINE ***') 


% *** calculation of required parameters *** 

rho=.002377*(-.000031 *PA+(-.002*temp+1.118)); 

% Hardwire rho for high and hot conditions 
if PA==4000 & temp==95 
rho=.00I92 
else 
end 

% *** first guess at rotor profile drag ( H force ) *** 

ifVinf< 16.9, 

Drotor=0; 

else 

Drotor=Vinf*(rho/.002377); 

end 


q=0.5*rho*Vinf''2; 

Adisk=pi*R''2; 

Vtip=omega*R; 

Dfuse=q*Afh; 

CDwing=CDowing+(CLwing''2/(ewing*pi*(bwing^2/Swing))); 

CDhoriz=CDohoriz+(CLhoriz''2/(.8*pi*(bhoriz''2/Shoriz))); 

CD vert=CDovert+(CLvert''2/(. 8 *pi * (bverU2/S vert))); 

Dwing=q*CDwing*Swing; 

Dhoriz=q*CDhoriz*Shoriz; 

Dvert=q*CDvert*Svert; 

Lwing=q*CLwing*Swing; 

Lhoriz=q*CLhoriz*Shoriz; 

Lvert=q*CLvert*Svert; 

LftotaI=Lwing+Lhoriz; 

DftotaI=Dfuse+Dwing+Dhoriz+Dvert; 

% *** thrust calculation *** 

% Check percentage of wing lift versus rotor thrust 

if heIochoice==2 & Vinf > 16.9 
cic 

LoverT=LftotaPGW; 

dispCThe wing is carrying '),disp(LoverT*IOO),disp('percent of the aircraft weight.') 

flag=I; 

while flag > 0 

dispC) 

dispCDo you want to change the percentage of aircraft weight carried by the wing? ') 
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disp('l. continue with current wing settings or 2. change percentage of aircraft weight 
carried by the wing by varying wing area') 

wingchoice=input('Enter a 1 or 2...'); 
while isempty(wingchoice), 
dispC ') 

dispCYou must enter a 1 or 2') 

dispCDo you want to change the percentage of aircraft weight carried by the wing? ') 
disp('l. continue with current wing settings or 2. change percentage of aircraft weight 
carried by the wing by varying wing area') 

wingchoice=input('Enter a 1 or 2...'); 
end 

if wingchoice~= 1 & wingchoice~=2 
dispC ') 

dispC *** Enter a 1 or 2 ***') 
elseif wingchoice==2 

wpercent=input('Enter desired percentage of aircraft weight carried by wing: '); 
lifterror=wpercent/100 * GW-Lftotal; 
if lifterror > 0 

while lifterror > .05*wpercent/100*GW 
Swing=Swing + 5; 

Lwing=q*CLwing*Swing; 

Lftotal=Lwing+Lhoriz; 
lifterror=wpercent/100 * GW-Lftotal; 
end 

lifterror=0; 

dispCNew wing area is '),disp(Swing),dispC fG2') 

dispCPress any key to continue ...') 

pause 

elseif lifterror < 0 

while lifterror < .05*wpercent/100*GW 
Swing=Swing - 5; 

Lwing=q*CLwing*Swing; 

Lftotal=LwingH-Lhoriz; 

lifterror=wpercent/100*GW-Lftotal; 

end 

lifterror=0; 

dispCNew wing area is '),disp(Swing),dispC fG2') 
dispCPress any key to continue ...') 
pause 
end 
flag=0; 

elseif wingchoice==l 

dispCWing settings unchanged') 
dispCPress any key to continue ...') 
pause 
flag=0; 
end 
end 


end 

% if required, define alphaT fixed quantity 
if alphaTfix==l 

alphaT=alphaTinput; 

else 

alphaT=atan2((DftotalH-Drotor), (GW-Lftotal)); 
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end 


% Determine thrust parameters for low speed flight 
ifVinf< 16.9, 

% Determine augmented vertical thrust capability 
clc 

dispCJANRAD is about to perform low speed performance analysis.') 
disp(") 
flag=l; 
while flag > 0 
disp(") 

dispCAre any auxiliary vertical thrust devices installed on the aircraft,') 

dispCi.e., Lift Fans or Direct Jet thrusters?') 

dispC) 

vertauxchoice=input('l. yes or 2. no'); 
while isempty(vertauxchoice), 
dispC ') 

dispCYou must enter a 1 or 2') 

dispCAre there any auxiliary vertical thrust devices installed on the aircraft,') 

dispCi.e., Lift Fans or Direct Jet thrusters?') 

dispC) 

vertauxchoice=inputCL yes or 2. no'); 
end 

if vertauxchoice~=l & vertauxchoice~=2 
dispC ') 

dispC *** Enter a 1 or 2 ***') 
elseif vertauxchoice== 1 

vertaux=inputCEnter amount of thrust in pounds provided by the vertical auxiliary thrust 

devices: '); 

GW=GW-vertaux; 

dispCInput parameters adjusted to reflect auxiliary vertical thrust.') 
flag=0; 

elseif vertauxchoice==2 
vertaux=0; 

dispCNo vertical thrust devices installed') 
dispCPress any key to continue ...') 
pause 
flag=0; 
end 
end 

T=(l+(0.3*Afv/Adisk))*GW; 

CT=T/(Adisk*rho*Vtip''2); 

% Determine thrust parmeters for nominal flight speeds 
else 

T=(GW-Lftotal)/cos(alphaT); 

CT=T/(Adisk*rho*Vtip''2); 

vertaux=0; 


% Added for Compound/RVR auxiliary thrust analysis 
if helochoice==2 & Vinf > 16.9 
clc 

if T auxchoice== 1 

T aux=(Dfuse+Dwing+Dhoriz+Dvert); 

% if required, define alphaT fixed quantity 
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if alphaTfix==l 

alphaT=alphaTinput; 

else 

alphaT=atan2((Dftotal+Drotor), (GW-Lftotal)); 
end 
else 

if Taux<((Dfuse+Dwing+Dhoriz+Dvert)-sin(alphaT)*T) 

dispC*** Auxiliary thrust will not overcome drag at this airspeed, ***') 
dispC Total drag = '),disp(Dfuse+Dwing+Dhoriz+Dvert) 

dispC Thrust deficit ='),disp(Dfuse+Dwing+Dhoriz+Dvert-sin(alphaT)*T-Taux) 
flag=l; 
while flag > 0 
dispC) 

dispCDo you want to change the auxiliary thrust value? ') 

disp(T. continue with current auxiliary thrust value or 2. set auxiliary thrust to equal 

drag') 

Tauxchoice=input('Enter a 1 or 2...'); 
while isempty(Tauxchoice), 
dispC ') 

dispC You must enter a 1 or 2') 
dispC) 

dispCDo you want to change the auxiliary thrust value? ') 
dispC) 

dispC 1. continue with current auxiliary thrust value or 2. set auxiliary thrust to equal 

drag') 

Tauxchoice=input('Enter a 1 or 2...'); 
end 

if Tauxchoice~=l & Tauxchoice~=2 
disp('') 

disp(' *** Enter a 1 or 2 ***') 
elseif Tauxchoice==2 

T aux=(Dfuse+Dwing+Dhoriz+Dvert); 
flag=0; 

disp('Auxiliary thrust change to equal total drag') 
disp(") 

disp('Press any key to continue ...') 
pause 

elseif T auxchoice== 1 

disp('Auxiliary thrust deficit is '),disp(Dfuse+Dwing+Dhoriz+Dvert-sin(alphaT)*T- 

Taux),disp(' lbs.') 

disp('*** Results may be erroneous. ***') 
disp(") 

disp('Press any key to continue ...') 
pause 
flag=0; 
end 
end 
end 

% define alphaT fixed quantity 
if alphaTfix==l 

alphaT=alphaTinput; 

else 

alphaT=atan2((Dftotal+Drotor), (GW-Lftotal)); 
end 
end 
end 


114 



Dftotal=(Dfuse+Dwing+Dhoriz+Dvert)-Taux; 
mu=Vinf*cos(alphaT)/V tip; 

optmu=.0107*(Vinf/1.69)-1.2; 


% if desired, iterate omega to obtain desired mu value 
if helochoice==2 
ifrvrCP==l 
ifmu< 1.0 
clc 

dispCThe advance ratio, mu, is'),disp(mu) 
disp(") 
flag=l; 
while flag > 0 

dispCThis value is too low for RVR operation, do you want JANRAD to optimize mu so 

the') 

dispCrotor system will operate in RVR mode?') 

disp('l. continue with current mu value or 2. change mu by automatically varying 
rotational speed'); 

muchoice=input('Enter a 1 or 2...'); 
while isempty(muchoice), 
dispC ') 

dispCYou must enter a 1 or 2') 

dispCDo you want JANRAD to optimize mu so the rotor system will operate in RVR 

mode? ') 

disp('l. continue with current mu value or 2. change mu by automatically varying 

rotational speed'); 

muchoice=input('Enter a 1 or 2...'); 
dispCOptimum mu value currently set to'),disp(optmu) 
end 

if muchoice~=l & muchoice~=2 
dispC ') 

dispC *** Enter a 1 or 2 ***') 
elseif muchoice==2 
muerror=mu-optmu; 
if muerror < 0 

while abs(muerror) > .01 
omega=omega - .001; 

Vtip=omega*R; 
mu=Vinf*cos(alphaT)A^ tip; 
muerror=mu-optmu; 
end 

muerror=0; 

dispCNew rotational velocity is '),disp(omega) 
dispC) 

dispCNew mu value is '),disp(mu) 
dispC) 

dispCPress any key to continue ...') 
pause 
else 
end 
flag=0; 

elseif muchoice==l 

dispCRotational velocity and mu value unchanged.') 
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disp(") 

dispCPress any key to continue ...') 
pause 
flag=0; 
end 
end 
else 
end 
else 
end 
else 
end 

ifmu<= 1.0&rvrCP==l 
clc 

dispCThe advance ratio is less than 1.0 and you have chosen RVR based') 
dispCcyclic pitch. JANRAD cannot complete computations using this input') 
dispCThe program will revert to conventional Higher Harmonic Feathering.') 
dispCPress any key to continue ...') 
pause 
rvrCP=2; 
else 
end 


if rvrCP==l & twoPinput==0 
clc 

dispCYou have selected RVR cyclic pitch, 2/rev cyclic input cannot be zero.') 
twoPinput=inputCPlease enter a non-zero 2/rev cyclic input in degrees or press <CONTROL H- 
C> to quit JANRAD'); 

while isempty(twoPinput), 
dispC') 

dispCYou must enter a value') 
twoPinput=inputC2/rev cyclic input (degs): '); 
while twoPinput==0 

dispCThe value cannot be zero either!') 

twoPinput=inputCPlease enter a non-zero 2/rev cyclic input or press <CONTROL + C> to 

quit JANRAD'); 

end 

end 

else 

end 

clc 

dispC *** ANALYSIS IN PROGRESS ***') 

% *** setup blade radius elements, azimuth elements, 

% induced velocity distributions, and determination 
% of coning angle and tip loss parameters *** 

B=l-(sqrt(2*CT)/b); 

Reff=B*R; 

Rbar=Reff-e; 
dr=(Reff-grip)/nbe; 
r=grip:dr;Reff-dr;,r=rH-dr/2; 
rTl=.7;,% *** first guess at rT *** 

RbarT=rTl*Rbar; 
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mblade=wblade/32.17; 

betao=asin((T/b*RbarT-(.5*(R-e)+e)*wblade)/((.5*(R-e)+e)^2*omega^2*mblade)); 

betat=twist*(0.7-(r/R)); 

psi=0:360/naz:360-360/naz;,psi=psi757.3; 

% *** induced velocity distribution *** 

ifVinf< 16.9, 

A=4*pi; 

Bv=(b/2)*omega*a*cblade; 

Tv=0; 

delT=T-Tv; 

while abs(delT) > .01*T 

thetav=twist*(0.7-(r/R))+thetao; 

C=(-b/2)*cblade*omega^2*r*a.*thetav; 

vi=(-Bv+sqrt(Bv^2-(4*A*C)))/(2*A); 

dTv=(b/2)*rho*((omega*r).^2)*a.*(thetav-(vi./(omega*r)))*cblade*dr; 

Tv=sum(dTv); 

delT=T-Tv; 

ifdelT<0, 

thetao=thetao-0.5*thetao*abs(delT/T); 

else 

thetao=thetao+0.5*thetao*abs(delT/T); 

end 

end 

else 

lamdaT=[l-2*mu*sin(alphaT) mu^2*(sin(alphaT)^2+l)-2*mu^3*sin(alphaT) 

muM*sin(alphaT)^2-CT^2/4]; 

lamdaT=max(real(roots(lamdaT))); 
vi=lamdaT* Vtip-V inf* sin(alphaT); 
vi=vi*ones(size(r)); 
end 

% *** first guess at theta *** 

thetalc=0.035*((0.0006e-3*Vinf''2+0.244e-3*Vinf)/0.105); 

thetals=-0.087*((0.0006e-3*Vinf''2+0.244e-3*Vinf)/0.105); 

theta2c=twoPinput/57.3; 

theta3s=threePinput/57.3; 

theta3Poffset=threePoffset/57.3; 

theta=thetao+theta 1 c. *cos(psi)+theta 1 s. * sin(psi); 

theta2Pr=theta2c.*cos(2*psi); 

theta3Pr=theta3s.*sin(3*psi+theta3Poffset); 

thetaXP=thetao+thetalc.*cos(psi)+thetals.*sin(psi)+theta2c.*cos(2*psi)+theta3s.*sin(3*psi+theta 

3Poffset); 

if rvrCP== 1 

thetaXP=thetao+thetao*sin(psi)+thetalc.*cos(psi)+thetals.*sin(psi)+theta2c.*cos(2*psi); 

thetaXP=thetao+thetao*sin(psi)-thetalc.*cos(psi)-thetals.*sin(psi)+theta2c.*cos(2*psi); 

else 

end 

% *** rotor trimming routine *** 

dispC ') 
dispC ') 
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dispC *** TRIMMING COLLECTIVE ***') 
k=I; 


ifrvrCP==I 

devTC=.20; 

devTC=.I; 

else 

devTC=.02; 

end 

errorO=(T*devTC)+1; 
while abs(errorO) > T*devTC; 

Tpsi=zeros(size(psi)); 

thrcalc 

errorO=T-(mean(Tpsi) *b); 
if errorO < -T*devTC 

thetao=thetao-0.35*thetao*abs(L5*error0/T)*abs(I-mu);% added abs for RVR operation 
elseif errorO > T*devTC 

thetao=thetao+0.35*thetao*abs(L5*error0/T)*abs(I-mu);% added abs for RVR operation 
end 

theta=thetao+theta I c. *cos(psi)+theta I s. * sin(psi); 
theta2Pr=theta2c.*cos(2*psi); 
theta3Pr=theta3 s. * sin(3 *psi+theta3Poffset); 

thetaXP=thetao+thetaIc.*cos(psi)+thetaIs.*sin(psi)+theta2c.*cos(2*psi)+theta3s.*sin(3*psi+theta3Poffset) 


if rvrCP== I 

thetaXP=thetao+thetao*sin(psi)+thetaIc.*cos(psi)+thetaIs.*sin(psi)+theta2c.*cos(2*psi); 

thetaXP=thetao+thetao*sin(psi)-thetaIc.*cos(psi)-thetaIs.*sin(psi)+theta2c.*cos(2*psi); 

else 

end 

if k > I, 

if abs(errorO) >= abs(errorl), 
cic 

dispC ') 
dispC ') 

dispCThis configuration will not trim') 
dispCTry a lower airspeed or a new design') 
dispC ') 

dispCProgram execution terminated') 
dispC ') 

errorC*** END OE PROGRAM ***') 
end 
end 

errorl=error0; 
ifrvrCP==I 
errorO=(T*devTC)-I; 
end 
end 

dispC ') 

dispC *** TRIMMING CYCLIC ***') 
t0=clock; 
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k=l; 

ifrvrCP==2 

errorO=(((T/b)*rTl *(R-grip))*.04)+1; 
while errorO > ((T/b)*rTl*(R-grip))*.04 
time=etime(clock,tO); 

if time > 15, 
dispC ') 

dispCstill trimming') 
tO=clock; 
end 

Mpsi(:,k)=zeros(size(psi)); 

tmcalc 

theta= [theta theta(:,k)]; 
theta2Pr= [theta2Pr theta2Pr(: ,k)]; 
the ta3 Pr= [ thetaS Pr thet a3 Pr(: ,k) ]; 
thetaXP=[thetaXP thetaXP(:,k)]; 

Mpsi=[Mpsi Mpsi(:,k)]; 

% *** calculation of initial dthetadM *** 

ifk<2, 

theta(:,k+l)=theta(:,k)+0.25/57.3; 

theta2Pr(:,k+l)=theta2Pr(;,k)+0.25/57.3; 

thetaXP(:,k+l)=thetaXP(:,k)+0.25/57.3; 

Mpsi(: ,k+1 )=zeros(size(psi)); 

k=k+l; 

tmcalc 

k=k-l; 

dthetadM=(thetaXP(: ,k+1 )-thetaXP(:,k))./(Mpsi(: ,k+1 )-Mpsi(: ,k)); 
end 

% *** calculation of M first harmonic parameters *** 

M1 c=2* sum(Mpsi(: ,k). *cos(psi))/naz; 

M1 s=2* sum(Mpsi(: ,k). * sin(psi))/naz; 

% *** removal of first harmonic terms from Mpsi *** 

Mpsi(:,k+l)=Mpsi(:,k)-Mlc.*cos(psi)-Mls.*sin(psi); 
delM=Mpsi(: ,k+1 )-Mpsi(:,k); 
errorO=max(delM)-min(delM); 
ifk> 1, 

if errorO > error 1, 
clc 

dispC ') 
dispC ') 

dispCThis configuration will not trim') 
dispCTry a lower airspeed or a new design') 
dispC ') 

dispCProgram execution terminated') 
dispC ') 

errorC*** END OF PROGRAM ***') 
end 
end 
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error l=eirorO; 


% *** calculation of new theta *** 

delM=0.5*(l-mu)*delM; 

theta(: ,k+1 )=theta(: ,k)+(dthetadM. *delM); 

theta2Pr(; ,k+1 )=theta2Pr(: ,k)+(dthetadM. *delM); 

theta3Pr(; ,k+1 )=theta3Pr(: ,k)+(dthetadM. *delM); 

thetaXP(:,k+l)=thetaXP(;,k)+(dthetadM.*delM); 

if errorO <= ((T/b)*rTl*(R-grip))*.04, 
theta 1 c=2* sum(theta(: ,k). *cos(psi))/naz; 
theta 1 s=2* sum(theta(: ,k). * sin(psi))/naz; 

% *** 2P coefficient *** 

theta2c=sum(ones(size(theta(:,k))).*twoPinput/57.3)/naz; 

% *** 3P coeeficient *** 

theta3s=sum(ones(size(theta(:,k))).*threePinput/57.3)/naz; 

else 

theta 1 c=2* sum(theta(: ,k+1). *cos(psi))/naz; 
theta 1 s=2* sum(theta(: ,k+1). * sin(psi))/naz; 

% *** 2P coefficient *** 

theta2c=sum(ones(size(theta(:,k+l))).*twoPinput/57.3)/naz; 

% *** 3P coefficient *** 

theta3s=sum(ones(size(theta(:,k+l))).*threePinput/57.3)/naz; 

end 

% *** calculate theta with IP terms only *** 

theta(: ,k+1 )=thetao+theta 1 c. *cos(psi)+theta 1 s. * sin(psi); 

% *** add 2P and 3P terms, result thetaXP *** 
theta2Pr(; ,k+1 )=theta2c. *cos(2*psi); 
theta3Pr(; ,k+1 )=theta3 s. * sin(3 *psi+theta3Poffset); 

thetaXP(:,k+l)=thetao+thetalc.*cos(psi)+thetals.*sin(psi)+theta2c.*cos(2*psi)+theta3s.*sin(3*psi+theta3 

Poffset); 

ifrvrCP==l 


thetaXP(:,k+l)=thetao+thetao*sin(psi)+thetalc.*cos(psi)+thetals.*sin(psi)+theta2c.*cos(2*psi); 

thetaXP=thetao+thetao*sin(psi)-thetalc.*cos(psi)-thetals.*sin(psi)+theta2c.*cos(2*psi); 

else 

end 

% *** calculation of new dthetadM *** 
theta=[theta theta(:,k+l)]; 
theta2Pr= [theta2Pr theta2Pr(: ,k+1)]; 
theta3Pr=[theta3Pr theta3Pr(:,k+l)]; 
thetaXP= [thetaXP thetaXP(:,k+l)]; 

Mpsi=[Mpsi Mpsi(:,k+1)]; 

theta(: ,k+2)=theta(: ,k)+0.25/57.3; 
theta2Pr(; ,k+2)=theta2Pr(: ,k)+0.25/57.3; 
theta3Pr(; ,k+2)=theta3Pr(: ,k)+0.25/57.3; 
thetaXP(:,k+2)=thetaXP(;,k)+0.25/57.3; 

Mpsi(: ,k+2)=zeros(size(Mpsi(: ,k+1))); 
k=k+2; 
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tmcalc 

k=k-2; 

dthetadM=(thetaXP(:,k+2)-thetaXP(:,k))./(Mpsi(:,k+2)-Mpsi(:,k)); 

k=k+l; 

end 

else 

Mpsi(:,k)=zeros(size(psi)); 

tmcalc 

theta=[theta theta(:,k)]; 
theta2Pr=[theta2Pr theta2Pr(:,k)]; 
theta3Pr=[theta3Pr theta3Pr(:,k)]; 
thetaXP=[thetaXP thetaXP(;,k)]; 

Mpsi=[Mpsi Mpsi(:,k)]; 

% *** calculation of initial dthetadM *** 

ifk<2, 

theta(:,k+1 )=theta(: ,k)+0.25/57.3; 

theta2Pr(:,k+l)=theta2Pr(;,k)+0.25/57.3; 

thetaXP(:,k+l)=thetaXP(:,k)+0.25/57.3; 

Mpsi(:,k+l)=zeros(size(psi)); 

k=k+l; 

tmcalc 

k=k-l; 

dthetadM=(thetaXP(: ,k+1 )-thetaXP(:,k))./(Mpsi(: ,k+1 )-Mpsi(: ,k)); 

end 

end 

dispC ') 

dispC *** ADJUSTING COLLECTIVE ***') 
dispC ') 

theta=theta(:,k); 
theta2Pr=theta2Pr(; ,k); 
theta3Pr=theta3Pr(; ,k); 
thetaXP=thetaXP(;,k); 

k=l; 


ifrvrCP==l 

devAC=.l; 

devAC=.03; 

else 

devAC=.01; 

end 

errorO=(T*devAC)+1; 
while abs(errorO) > T*devAC 
Tpsi=zeros(size(psi)); 
thrcalc 

errorO=T-(mean(Tpsi) *b); 
if errorO < -T*devAC, 

thetao=thetao-0.25*thetao*abs(L25*error0/T)*abs(l-mu); 
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elseif errorO > T*devAC, 

thetao=thetao+.025*thetao*abs(1.25*error0/T)*abs(l-mu); 

end 

theta=thetao+theta 1 c. *cos(psi)+theta Is.* sin(psi); 
theta2Pr=theta2c.*cos(2*psi); 
theta3Pr=theta3 s. * sin(3 *psi+threePoffset); 

thetaXP=thetao+thetalc.*cos(psi)+thetals.*sin(psi)+theta2c.*cos(2*psi)+theta3s.*sin(3*psi+theta3Poffset) 


if rvrCP== 1 

thetaXP=thetao+thetao*sin(psi)+thetalc.*cos(psi)+thetals.*sin(psi)+theta2c.*cos(2*psi); 

thetaXP=thetao+thetao*sin(psi)-thetalc.*cos(psi)-thetals.*sin(psi)+theta2c.*cos(2*psi); 

else 

end 

if k > 1, 

if abs(errorO) > abs(errorl), 
clc 

dispC ') 
dispC ') 

dispCThis configuration will not trim') 
dispCTry a lower airspeed or a new design') 
dispC ') 

dispCProgram execution terminated') 
dispC ') 

errorCEND OF PROGRAM') 
end 
end 

errorl=errorO; 
k=k+l; 
ifrvrCP==l 
errorO=(T*devAC)-1; 
end 
end 

% *** calculating drag moments *** 

DMpsi=zeros(size(psi)); 
dmcalc 

% *** calculating rotor H force *** 

ifVinf< 16.9, 

Hrotor=0; 
dT=[dT ddT]; 
dN=[dN ddN]; 
dD=[dD ddD]; 
else 

dT=[dT ddT]; 
dN=[dN ddN]; 
dD=[dD ddD]; 
for i=l:length(r)+l, 

Hlc(i)=2*sum(dT(:,i).*cos(psi))/naz; 

Hls(i)=2*sum(dT(:,i).*sin(psi))/naz; 

end 

if rvrCP==2 

Hrotor=(((b*cos(alphaT)/2)*(sum(Hls)-sin(betao)*sum(Hlc)))+Drotor)/2; 
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else 

Hrotor=0; 

end 

end 

% *** calculating new rT *** 

rT2=(((mean(Mpsi(:,length(Mpsi(l,:))-l))/mean(Tpsi))/R)+rTl)/2; 
% *** check rotor drag and rT, retrim rotor if needed *** 


while abs(Drotor-Hrotor) > 0.2*Hrotor I abs(rTl-rT2) > 0.5*rTl % increase tolerance to 50% 
if abs(Drotor-Hrotor) > 0.2*Hrotor, 
dispC') 

dispC * * * ADJUSTING ROTOR DRAG * * *') 

end 


ifrvrCP==l 

Drotor=0; 

else 

Drotor=Hrotor; 

end 


ifrvrCP==l 

devMT=.25; 

devMT=.045; 

else 

devMT=.015; 

end 

if abs(rTl-rT2) > devMT*rTl, 
dispC') 

dispC *** ADJUSTING MEAN THRUST LOCATION ***') 
end 

dispC ') 

dispC *** RETRIMMING ROTOR ***') 
dispC ') 

dT=dT(:,l:nbe); 

dN=dN(:,l:nbe); 

dD=dD(:,l:nbe); 

% *** recealculating parameters *** 

% if required, define alphaT fixed quantity 
if alphaTfix==l 

alphaT=alphaTinput; 

else 

alphaT=atan((Dftotal+Drotor)/(GW-Lftotal)); 

end 

mu=Vinf*cos(alphaT)A^ tip; 

ifVinf>= 16.9, 

T=(GW-Lftotal)/cos(alphaT); 
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CT=T/(Adisk*rho*Vtip''2); 

lamdaT=[l -2*mu*sin(alphaT) mu''2*(sin(alphaT)^2+l)-2*mu''3*sin(alphaT) 

muM*sin(alphaT)^2-CT^2/4]; 

lamdaT=max(real(roots(lamdaT))); 

vi=lamdaT*Vtip-Vinf*sin(alphaT); 

vi=vi*ones(size(r)); 

end 

B=l-(sqrt(2*CT)/b); 

Reff=B*R; 

Rbar=Reff-e; 

dr=(Reff-grip)/nbe; 

r=grip:dr:Reff-dr;,r=r+dr/2; 

RbarT=rT2*Rbar; 

betao=asin((T/b*RbarT-(.5*(R-e)+e)*wblade)/((.5*(R-e)+e)''2*omega''2*mblade)); 

% *** trimming collective *** 

tO=clock; 

k=l; 

errorO=(T*devTC)+l; 

while abs(errorO) > T*devTC 
Tpsi=zeros(size(psi)); 
thrcalc 

errorO=T-(mean(Tpsi)*b); 
if errorO < -T*devTC, 

thetao=thetao-0.35*thetao*abs(1.5*error0/T)*abs(l-mu); 
elseif errorO > T*devTC, 

thetao=thetao+0.35*thetao*abs(1.5*error0/T)*abs(l-mu); 

end 

theta=thetao+theta 1 c. *cos(psi)+theta Is.* sin(psi); 
theta2Pr=theta2c.*cos(2*psi); 
theta3Pr=theta3 s. * sin(3 *psi+threePoffset); 

thetaXP=thetao+thetalc.*cos(psi)+thetals.*sin(psi)+theta2c.*cos(2*psi)+theta3s.*sin(3*psi+theta3Poffset) 


if rvrCP== 1 

thetaXP=thetao+thetao*sin(psi)+thetalc.*cos(psi)+thetals.*sin(psi)+theta2c.*cos(2*psi); 

thetaXP=thetao+thetao*sin(psi)-thetalc.*cos(psi)-thetals.*sin(psi)+theta2c.*cos(2*psi); 

else 

end 

if k > 1, 

if abs(errorO) > abs(errorl), 
clc 

dispC ') 
dispC ') 

dispCThis configuration will not trim') 
dispCTry a lower airspeed or a new design') 
dispC ') 

dispCProgram execution terminated at line 746') 
dispC ') 

error)'*** END OF PROGRAM ***') 
end 
end 

error l=errorO; 
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k=k+l; 
ifrvrCP==l 
errorO=(T*devTC)-1; 
end 
end 

% *** trimming cyclic *** 


ifrvrCP==l 

devCY=.2; 

devCY=.12; 

else 

devCY=.04; 

end 

k=l; 

errorO=(((T/b)*rT2*(R-grip))*devCY)+l; 


if rvrCP==2 

while errorO > ((T/b)*rT2*(R-grip))*devCY 
time=etime(clock,tO); 
if time > 15, 

dispCstill trimming...') 
dispC ') 
tO=clock; 
end 

Mpsi(:,k)=zeros(size(psi)); 

tmcalc 

theta= [theta theta(:,k)]; 
theta2Pr= [theta2Pr theta2Pr(: ,k)]; 
theta3Pr=[theta3Pr theta3Pr(:,k)]; 
thetaXP=[thetaXP thetaXP(:,k)]; 

Mpsi=[Mpsi Mpsi(:,k)]; 

% *** calculation of initial dthetadM *** 

ifk<2, 

theta(: ,k+1 )=theta(: ,k)+0.25/57.3; 

theta2Pr(:,k+l)=theta2Pr(:,k)+0.25/57.3; 

theta3Pr(:,k+l)=theta3Pr(:,k)+0.25/57.3; 

thetaXP(:,k+l)=thetaXP(:,k)+0.25/57.3; 

Mpsi(:,k+l)=zeros(size(psi)); 

k=k+l; 

tmcalc 

k=k-l; 

dthetadM=(thetaXP(:,k+l)-thetaXP(:,k))./(Mpsi(:,k+l)-Mpsi(:,k)); 

end 

% *** calculation of M first harmonic parameters *** 

Mlc=2*sum(Mpsi(:,k).*cos(psi))/naz; 

M1 s=2* sum(Mpsi(: ,k). * sin(psi))/naz; 

% *** removal of first harmonic terms from Mpsi *** 

Mpsi(:,k+l)=Mpsi(:,k)-Mlc.*cos(psi)-Mls.*sin(psi); 
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delM=Mpsi(:,k+1 )-Mpsi(: ,k); 
errorO=max(delM)-min(delM); 
ifk> 1, 

if errorO > errorl, 
clc 

dispC ') 

dispCThis configuration will not trim') 
dispCTry a lower airspeed or a new design') 
dispC ') 

dispCProgram Execution Terminated at line 814') 
dispC ') 

error)'*** END OE PROGRAM ***') 
end 
end 

error l=errorO; 

% *** calculation of new theta *** 

delM=0.5*(l-mu)*delM; 

theta): ,k+1 )=theta): ,k)+)dthetadM. *delM); 

theta2Pr):,k+l)=theta2Pr);,k)+)dthetadM.*delM); 

theta3Pr):,k+l)=theta3Pr);,k)+)dthetadM.*delM); 

thetaXP): ,k+1 )=thetaXP): ,k)+)dthetadM.*delM); 

if errorO <= ))T/b)*rT2*)R-grip))*.04, 
theta 1 c=2* sum)theta): ,k). *cos)psi))/naz; 
theta 1 s=2* sum)theta): ,k). * sin)psi))/naz; 

% 2P/3P coefficients 

theta2c=sum)ones)size)theta):,k))).*twoPinput/57.3)/naz; 

theta3s=sum)ones)size)theta):,k))).*threePinput/57.3)/naz; 

else 

thetalc=2*sum)theta):,k+l).*cos)psi))/naz; 
theta 1 s=2* sum)theta): ,k+1). * sin)psi))/naz; 

% 2P/3P coefficients 

theta2c=sum)ones)size)theta):,k+l))).*twoPinput/57.3)/naz; 
theta3 s=sum)ones)size)theta): ,k+1))). * threePinput/57.3)/naz; 


% calculate theta with IP terms only 

theta): ,k+1 )=thetao+theta 1 c. *cos)psi)+theta 1 s. * sin)psi); 

% add 2P and 3P terms, unit value only for theta2c, result thetaXP 
theta2Pr): ,k+1 )=theta2c. *cos)2*psi); 
theta3Pr):,k+l)=theta3s.*sin)3*psi+threePoffset); 

thetaXP):,k+l)=thetao+thetalc.*cos)psi)+thetals.*sin)psi)+theta2c.*cos)2*psi)+theta3s.*sin)3*psi+theta3 
Poffset); 


if rvrCP==l 


thetaXP):,k+l)=thetao+thetao*sin)psi)+thetalc.*cos)psi)+thetals.*sin)psi)+theta2c.*cos)2*psi); 

thetaXP):,k+l)=thetao+thetao*sin)psi)-thetalc.*cos)psi)- 
theta Is.* sin)psi)+theta2c .*cos)2*psi); 
else 
end 

% *** calculation of new dthetadM *** 
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theta=[theta theta(:,k+l)]; 
theta2Pr= [theta2Pr theta2Pr(: ,k+1)]; 
theta3Pr=[theta3Pr theta3Pr(:,k+1)]; 
thetaXP=[thetaXP thetaXP(:,k+l)]; 

Mpsi=[Mpsi Mpsi(:,k+1)]; 
theta(; ,k+2)=zeros(size(Mpsi(; ,k+1))); 
theta2Pr(: ,k+2)=zeros(size(Mpsi(: ,k+1))); 
theta3Pr(: ,k+2)=zeros(size(Mpsi(: ,k+1))); 
thetaXP(: ,k+2)=zeros(size(Mpsi(: ,k+1))); 

k=k+2; 

tmcalc 

k=k-2; 

% replace IP dthetadM with 2P dthetadM 
% dthetadM=(theta(: ,k+2)-theta(: ,k)) ./(Mpsi(: ,k+2)-Mpsi(: ,k)); 

dthetadM=(thetaXP(:,k+2)-thetaXP(:,k))./(Mpsi(:,k+2)-Mpsi(:,k)); 

k=k+l; 

end 

else 

Mpsi(:,k)=zeros(size(psi)); 

tmcalc 

theta= [theta theta(:,k)]; 
theta2Pr= [theta2Pr theta2Pr(: ,k)]; 
theta3Pr=[theta3Pr theta3Pr(:,k)]; 
thetaXP=[thetaXP thetaXP(:,k)]; 

Mpsi=[Mpsi Mpsi(:,k)]; 

% *** calculation of initial dthetadM *** 

ifk<2, 

theta(:,k+l)=theta(:,k)+0.25/57.3; 

theta2Pr(:,k+l)=theta2Pr(:,k)+0.25/57.3; 

theta3Pr(:,k+l)=theta3Pr(:,k)+0.25/57.3; 

thetaXP(:,k+l)=thetaXP(:,k)+0.25/57.3; 

Mpsi(: ,k+1 )=zeros(size(psi)); 

k=k+l; 

tmcalc 

k=k-l; 

dthetadM=(thetaXP(: ,k+1 )-thetaXP(:,k))./(Mpsi(: ,k+1 )-Mpsi(: ,k)); 
end 
end 

% *** re trimming collective *** 

theta=theta(:,k); 
theta2Pr=theta2Pr(: ,k); 
theta3Pr=theta3Pr(: ,k); 
thetaXP=thetaXP(:,k); 
k=l; 

errorO=(T*devAC)+l; 
while abs(errorO) > T*devAC 
Tpsi=zeros(size(psi)); 
thrcalc 
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errorO=T-(mean(Tpsi)*b); 
if errorO < -T*devAC, 

thetao=thetao-0.25*thetao*abs(1.25*error0/T)*abs(l-mu); 
elseif errorO > T*devAC, 

thetao=thetao+0.25*thetao*abs(1.25*error0/T)*abs(l-mu); 

end 

theta=thetao+theta 1 c. *cos(psi)+theta Is.* sin(psi); 
theta2Pr=theta2c.*cos(2*psi); 
theta3Pr=theta3 s. * sin(3 *psi+threePoffset); 

thetaXP=thetao+thetalc.*cos(psi)+thetals.*sin(psi)+theta2c.*cos(2*psi)+theta3s.*sin(3*psi+theta3Poffset) 


if rvrCP== 1 

thetaXP=thetao+thetao*sin(psi)+thetalc.*cos(psi)+thetals.*sin(psi)+theta2c.*cos(2*psi); 

thetaXP=thetao+thetao*sin(psi)-thetalc.*cos(psi)-thetals.*sin(psi)+theta2c.*cos(2*psi); 

else 

end 

if k > 1, 

if abs(errorO) > abs(errorl), 
clc 

dispC ') 

dispCThis configuration will not trim') 
dispCTry a lower airspeed or a nex design') 
dispC ') 

error('END OF PROGRAM at line 933') 
end 
end 

error l=errorO; 
k=k+l; 
ifrvrCP==l 
errorO=(T*devAC)-1; 
end 
end 

% *** recalculating rotor H force *** 

ifVinf< 16.9, 

Hrotor=0; 
dT=[dT ddT]; 
dN=[dN ddN]; 
dD=[dD ddD]; 
else 

dT=[dT ddT]; 
dN=[dN ddN]; 
dD=[dD ddD]; 
for i= 1 :length(r)+1, 

Hlc(i)=2*sum(dT(:,i).*cos(psi))/naz; 

Hls(i)=2*sum(dD(:,i).*sin(psi))/naz; 

end 


ifrvrCP==l 

Hrotor=0; 

else 

Hrotor=(((b*cos(alphaT)/2)*(sum(Hls)-sin(betao)*sum(Hlc)))+Drotor)/2; 

end 
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end 


% * * * recalculating rT * * * 
rTl=rT2; 

rT2=(((mean(Mpsi(:,length(Mpsi(l,:))-l))/mean(Tpsi))/R)+rTl)/2; 

end 

% *** recalculating drag moments *** 

dT=dT(;,l;nbe); 

dN=dN(:,l;nbe); 

dD=dD(:,l:nbe); 

DMpsi=zeros(size(psi)); 

dmcalc 

dT=[dT ddT]; 

dN=[dN ddN]; 

dD=[dD ddD]; 

clc 

dispC ') 

dispC *** ROTOR TRIMMED ***') 

dispC press any key to continue...') 

pause 
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APPENDIX!. 0012CLCD.M 


function [CL,CD]=ool2clcd(alpha, mach) 

% 

% ool2clcd calculates CL and CD for the NACA 0012 
% airfoil given angle of attack in radians and the 
% local Mach number: 

% 

% [CL,CD]=ool2clcd(alpha, mach) 

% 

% Both 'alpha' and 'mach' are intended to be vectors 
% the elements of which correspond to the rotor blade 
% radial stations of interest in a blade element analysis. 
% All equations are based on Ray Prouty's treatment of 
% the 0012 in his text. 


CL=zeros(size(alpha)); 
CD=zeros(size(alpha)); 
a=alpha* 180/pi; 
mach=abs(mach); 
aL = 15 - 16.*mach; 
aD = 17 - 23.4.*mach; 

K1 = 0.0233 + 0.342.*(mach.^7.15); 
K2 = 2.05 - 0.95.*mach; 


% CL for Mach numbers < 0.725 and AOA inside +/- 20 deg: 

chk=(mach<0.725 & a>=0 & a<=aL); 

CL=CL+chk.*((0.L/sqrt(l-mach.''2) - 0.0L*mach).*a); 

chk=(mach<0.725 & a>aL & a<=20); 

CL=CL+chk.*((0.L/sqrt(l-mach.^2) - 0.0L*mach).*a - KL*(a-aL).''K2); 
chk=(mach<0.725 & a>=-20 & a<-aL); 

CL=CL-chk.*((0.L/sqrt(l-mach.^2) - 0.0L*mach).*abs(a) - KL*(abs(a)-aL).^K2); 

chk=(mach<0.725 & a>=-aL & a<0); 

CL=CL-chk.*((0.L/sqrt(l-mach.^2) - 0.0L*mach).*abs(a)); 


% CL for Mach numbers > 0.725 and AOA inside +/- 20 deg: 

chk=(mach>=0.725 & a>=0 & a<=aL); 

CL=CL+chk.*((0.677 - 0.744.*mach).*a); 

chk=(mach>=0.725 & a>aL & a<=20); 

CL=CL+chk.*((0.677 - 0.744.*mach).*a - (0.0575-0.144.*(mach-0.725).''0.44).*(a-aL).^(K2)); 

chk=(mach>=0.725 & a<0 & a>=-aL); 

CL=CL-chk.*((0.677 - 0.744.*mach).*abs(a)); 
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chk=(mach>=0.725 & a<-aL & a>=-20); 

CL=CL-chk.*((0.677 - 0.744.*mach).*abs(a) - (0.0575-0.144.*(mach-0.725).^0.44).*(abs(a)- 
aL).^K2)); 


% CL for all Mach numbers and AOA outside +!- 20deg: 

chk=(a>20&a<=161); 

CL=CLH-chk.*(1.15.*sin(2.*alpha)); 

chk=(a>161 & a<=173); 

CL=CLH-chk.*(-0.7); 

chk=(a>173 & a<=180); 

CL=CLH-chk.*(0.1.*(a-180)); 


chk=(a>=-180&a<=-173); 

CL=CLH-chk.*(0.1.*(aH-180)); 

chk=(a>-173 & a<=-161); 
CL=CLH-chk.*(0.7); 

chk=(a>-161 & a<-20); 
CL=CLH-chk.*(1.15.*sin(2.*alpha)); 


% CD for Mach numbers < 0.725 and AOA inside +!- 20 deg: 
chk=(mach<0.725 & a>=0 & a<=aD); 

CD=CDH-chk.*(0.0081 h- (-350.*a h- 396.*a.^2 - 63.3.*a.^3 h- 3.66.*a.M).*10.^(-6)); 
chk=(mach<0.725 & a>aD & a<=20); 

CD=CDH-chk.*((0.0081 h- (-350.*a h- 396.*a.^2 - 63.3.*a.^3 h- 3.66.*a.M).*10.^-6)) h- 
0.00066.*(a-aD).^2.54); 

chk=(mach<0.725 & a<0 & a>=-aD); 

CD=CDH-chk.*(0.0081 h- (-350.*abs(a) h- 396.*a.^2 - 63.3.*abs(a).^3 h- 3.66.*a.M).*10.^(-6)); 
chk=(mach<0.725 & a<-aD & a>=-20); 

CD=CDH-chk.*((0.0081 h- (-350.*abs(a) h- 396.*a.^2 - 63.3.*abs(a).^3 h- 3.66.*a.M).*10.^(-6)) h- 
0.00066.*(abs(a)-aD).^2.54); 


% CD for Mach numbers > 0.725 and AOA inside +!- 20 deg: 
chk=(mach>=0.725 & a>=0 & a<=20); 

CD=CDH-chk.*((0.0081 h- (-350.*a h- 396.*a.^2 - 63.3.*a.^3 h- 3.66.*a.M).*10.^-6)) h- 
0.00035.*a.''2.54 h- 21.*(mach-0.725).''3.2); 

chk=(mach>=0.725 & a<0 & a>=-20); 

CD=CDH-chk.*((0.0081 h- (-350.*abs(a) h- 396.*a.^2 - 63.3.*abs(a).^3 h- 3.66.*a.M).*10.^(-6)) h- 
0.00035.*abs(a).^2.54 H- 21.*(mach-0.725).^3.2); 
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% CD for all Mach numbers and AOA outside +/- 20deg: 


chk=(a>20 & a<=180); 

CD=CD+chk.*(1.03 - 1.02.*cos(2.*alpha)); 

chk=(a>=-180 & a<-20); 
CD=CD+chk.*(1.03 - 1.02.*cos(2.*alpha)); 
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APPENDIX J. HH02CLCD.M 


function [CL,CD]=hh02clcd(alpha) 

% hh02clcd calculates CL and CD for an HH-02 airfoil 
% given angle of attack (alpha) in radians and Mach number (Mach) 

% [CL,CD]=hh02clcd(alpha,mach) 

CL=zeros(size(alpha)); 

CD=zeros(size(alpha)); 
a=alpha* 180/pi; 

chkl=(a>=20&a<=180); 

CL=CL+chkl.*(0.42541+0.026863*a+5.5988e-4*a.^2-2.1493e-5*a.''3+1.5932e-7*a.M-3.4659e 

10*a.''5); 

CD=CD+chkl.*(-0.7179+0.061213*a-5.9861e-4*a.''2+7.3708e-6*a.^3-6.6605e-8*a.M+1.913e- 

10*a.^5); 


chkl=(a>=-180 & a<=-50); 

CL=CL+chkl.*(-4.6183-0.1923*a-3.5554e-3*a.^2-3.3273e-5*a.^3-1.4528e-7*a.M-2.3003e- 

10*a.^5); 

CD=CD+chkl.*(2.7093e-2-2.1309e-2*a+2.0335e-4*a.^2+3.47e-7*a.^3-3.0586e-8*a.M-1.2584e 

10*a.^5); 

chkl=(a>-50 & a<-20); 

CL=CL+chkl .*(-2.5519-0.22847*a-9.5667e-3*a.^2-1.705 le-4*a.^3-1.0909e-6*a.M); 

CD=CDH-chkl.*(2.7093e-2-2.1309e-2*aH-2.0335e-4*a.^2H-3.47e-7*a.^3-3.0586e-8*a.M-1.2584e 

10*a.^5); 


chkl=(a>=-20 & a<=-10); 

CL=CLH-chkl.*(-0.2H-0.089*aH-0.0034*a.^2); 

CD=CDH-chkl.*(2.7093e-2-2.1309e-2*aH-2.0335e-4*a.^2H-3.47e-7*a.^3-3.0586e-8*a.M-1.2584e 

10*a.''5); 

chkl=(a<20&a>-10); 

CL=CLH-chkl.*(5.8766e-2H-1.3131e-l*aH-2.4742e-3*a.^2-5.303e-4*a.^3-1.5818e-5*a.MH-1.28e- 

6*a.^5); 

chk2=a<-4; 

chk2=chk2.*chkl; 

CD=CDH-chk2.*(1.3786H-0.916*aH-0.21396*a.^2H-2.0371e-2*a.^3H-7.0076e-4*a.M); 
chk2=(a>=-4 & a<=7); 
chk2=chk2.*chkl; 

CD=CDH-chk2.*(9.732e-3H-3.2326e-4*aH-1.4392e-4*a.^2-8.5073e-5*a.^3H-1.1826e- 

6*a.MH-1.5271e-6*a.^5); 

chk2=a>7; 

chk2=chk2.*chkl; 

CD=CDH-chk2.*(1.842e-l-5.7532e-2*aH-5.8043e-3*a.^2-1.2803e-4*a.^3); 
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APPENDIX K. RVRCLCD.M 


function [CL,CD]=rvrclcd(alpha,mach) 

%RVRCLCD calculates CL and CD for a 12% thickness Fairchild 
% RVR airfoil given angle of attack (alpha) in radians and 
% mach number (mach). Reference Fairchild Report HC144R1070 
% figs D-28 and D-12 

% JANRAD ver lP2P3Pbeta 

CL=zeros(size(alpha)); 

CD=zeros(size(alpha)); 

Mach_CL=[0 .26 .4 .6 .8 .9]; 

aoa_CL=[-180 -150 -120 -90 -30 -24 -20 -16 -12 -8 -4 0 4 8 12 16 20 24 30 90 120 150 160 170 
180]; 


cl_tab=[ 

0 0 

0 

0 

0 

0; 

0 

.7 

.7 

.7 .7 

■7; 


0 

.4 

.4 

.4 .4 

■4; 


0 

0 

0 

0 0 

0; 


0 

-.75 

-.75 

-.75 

-.75 

-.75; 

0 

-.89 

-.89 

-.89 

-.89 

-.89; 

0 

-.9 

-.9 

-.9 -. 

9 -.' 

9; 

0 

-1.0 

-1.0 

-1.0 

-1.0 

-1.0; 

0 

-1.05 

-1.05 -1.05 -1.05 -1.05; 

0 

-.85 

-.85 

-.85 

-.85 

-.85; 

0 

-.5 

-.5 

-.5 -. 

5 -., 

5; 

0 

.25 

.25 

.25 

.25 

.25; 

0 

.25 

.25 

.25 

.25 

.25; 

0 

.7 

.7 

.7 .5 

.55; 

0 

.9 

.9 

.9 .7 

.65; 

0 

1.52 

1.52 

: 1.52 

.78 

1.0; 

0 

1.65 

1.65 

1.65 

.7 

.95; 

0 

1.5 

1.5 

1.5 

1.2 

1.2; 

0 

1.7 

1.7 

1.7 

1.7 

1.7; 

0 

0 

0 

0 0 

0; 


0 

-.55 

-.55 

-.55 

-.55 

-.55; 

0 

-1.09 

-1.09 -1.09 -1.09 -1.09; 

0 

-.8 

-.8 

-.8 -. 

65 - 

•6; 

0 

-1.1 

-1.1 

-1.1 

-.7 

-.5; 

0 

0 

0 

0 0 

0]: 


Mach_CD=[0 .26 

.4 .6.; 

8 .9]; 



aoa_CD= 

[-180 -150 -120 -90 -30 -24 

-20 -16 -12 -8 -4 0 4 8 12 16 20 24 30 90 120 150 160 170 

180]; 






cd_tab=[ 

0 0 

0 

0 

.02 

.15; 

0 

.65 

.65 

.65 

.65 

.65; 

0 

1.2 

1.2 

1.2 

1.2 

1.2; 

0 

1.65 

1.65 

1.65 

1.65 

1.65; 

0 

.6 

.6 

.6 .6 

•6; 


0 

.5 

.5 

.5 .5 

.5; 


0 

.42 

.42 

.42 

.42 

.42; 

0 

.4 

.4 

.4 .4 

■4; 


0 

.05 

.05 

.2 . 

2 .22; 
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0 

0 

0 

.15 

.15 

■2; 

0 

0 

0 

0 

.05 

.12; 

0 

0 

0 

0 

.1 

■1; 

0 

0 

0 

0 

.05 

.12; 

0 

0 

0 

.15 

.15 

■2; 

0 

.05 

.05 

.2 

.2 

.22; 

0 

.4 

.4 

.4 

.4 

■4; 

0 

.42 

.42 

.42 

.42 

.42; 

0 

.5 

.5 

.5 

.5 

.5; 

0 

.6 

.6 

.6 

.6 

.6; 

0 

1.65 

1.65 1.65 1. 

65 1. 

0 

1.2 

1.2 

1.2 

1.2 

1.2; 

0 

.65 

.65 

.65 

.65 

.65; 

0 

.45 

.45 

.45 

.45 

.45; 

0 

.15 

.15 

.17 

.2 

.21; 

0 

0 

0 

0 

.05 

.10]; 


rvra=alpha.* 180/pi; 
forj = l;length(rvra) 
ifrvra(i) < -180 
rvra(j) = rvra(j) + 360; 
elseif rvra(j) >180 
rvra(j) = rvra(j) - 360; 
end 
end 

mach = abs(mach); 

CL = interp2(Mach_CL,aoa_CL,cl_tab,mach,rvra); 
CD = interp2(Mach_CD,aoa_CD,cd_tab,mach,rvra); 
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APPENDIX L. SC1095R8CLCD.M 


function [CL,CD]=scl095r8clcd(alpha,mach) 

%scl095r8 calculates CL and CD for a Sikorsky SC1095R8 airfoil 
%given angle of attack (alpha) in radians and mach number (mach) 
% [CL,CD]=sc 1095r8clcd(alpha,mach) 

CL=zeros(size(alpha)); 

CD=zeros(size(alpha)); 


Mach_CL=[0 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0]; 

aoa_CL=[-180 -178 -176 -174 -172 -170 -168 -166 -164 -162 -160 -158 -90 -22 -20 -18 -16 -14 - 
12 - 10-8 ... 

-6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 90 158 160 162 164 166 168 170 172 174 176 178 

180]; 

cl_tab=[0. 0. 0. 0. 0. 0. 0. 0. 0.; 

.205 .205 .205 .205 .205 .205 .205 .205 .205; 


.41 

.41 

.41 

.41 

.41 

.41 

.41 

.41 

.41; 

.6 

.6 

.6 

.6 .6 

.6 

.6 

.6 

.6; 


.77 

.77 

.77 

.77 

.77 

.77 

.77 

.77 

.77; 

.82 

.82 

.82 

.82 

.82 

.82 

.82 

.82 

.82; 

.82 

.82 

.82 

.82 

.82 

.82 

.82 

.82 

.82; 

.8 

.8 

.8 

00 

00 

.8 

.8 

.8 

.8; 


.76 

.76 

.76 

.76 

.76 

.76 

.76 

.76 

.76; 

.705 

.705 .705 .705 .705 .705 

.705 

.705 .705; 

.65 

.65 

.65 

.65 

.65 

.65 

.65 

.65 

.65; 

.65 

.65 

.65 

.65 

.65 

.65 

.65 

.65 

.65; 


-.0627 -.0627 -.0627 -.0627 -.0627 -.0627 -.0627 -.0627 -.0627; 
-.98 -.98 -.98 -.914 -.934 -.926 -.875 -.838 -.822; 

-.975 -.975 -.96 -.910 -.93 -.920 -.856 -.81 -.79; 

-.969 -.969 -.962 -.906 -.926 -.914 -.838 -.782 -.758; 

-.963 -.963 -.966 -.902 -.922 -.908 -.819 -.754 -.726; 

-1.07 -1.07 -.824 -.803 -.805 -.88 -.8 -.726 -.694; 

-.718 -.718 -.532 -.528 -.66 -.83 -.79 -.698 -.662; 

-.366 -.366 -.24 -.4 -.61 -.78 -.81 -.67 -.630; 

-.245 -.245 -.3 -.33 -.55 -.74 -.75 -.666 -.622; 

-.39 -.39 -.44 -.32 -.52 -.68 -.69 -.663 -.615; 

-.4 -.4 -.42 -.44 -.47 -.58 -.47 -.486 -.428; 

-.185 -.185 -.185 -.196 -.199 -.255 -.25 -.31 -.24; 

.029 .029 .048 .048 .072 .07 .07 -.15 -.05; 

.244 .244 .282 .292 .343 .395 .35 .138 .2; 

.459 .459 .515 .536 .614 .72 .56 .39 .449; 

.673 .673 .749 .78 .84 .83 .705 .64 .7; 

.888 .888 .983 .96 .91 .882 .805 .765 .806; 

1.103 1.103 1.17 1.01 .946 .92 .842 .81 .85; 

1.25 1.25 1.13 .96 1. .924 .845 .829 .865; 

1.1 1.1 1.03 1.07 1.053 .928 .848 .848 .88; 

.98 .98 .96 1.06 1.075 .92 .860 .867 .895; 

.982 .982 .966 1.07 1.064 .9 .880 .886 .91; 

.984 .984 .972 1.065 1.053 .9 .900 .905 .925; 

.987 .987 .979 1.06 1.042 .92 .920 .924 .94; 

.0627 .0627 .0627 .0627 .0627 .0627 .0627 .0627 .0627; 
-.66 -.66 -.66 -.66 -.66 -.66 -.66 -.66 -. 66 ; 

-.655 -.655 -.655 -.655 -.655 -.655 -.655 -.655 -.655; 
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-.685 -.685 -.685 -.685 -.685 -.685 -.685 -.685 -.685; 
-.73 -.73 -.73 -.73 -.73 -.73 -.73 -.73 -.73; 

-.77 -.77 -.77 -.77 -.77 -.77 -.77 -.77 -.77; 


-.805 -.805 -.805 -.805 -.805 -.805 -.805 -.805 -.805; 
-.79 -.79 -.79 -.79 -.79 -.79 -.79 -.79 -.79; 

-.61 -.61 -.61 -.61 -.61 -.61 -.61 -.61 -.61; 

-.42 -.42 -.42 -.42 -.42 -.42 -.42 -.42 -.42; 

-.21 -.21 -.21 -.21 -.21 -.21 -.21 -.21 -. 21 ; 

0 . 0 . 0 . 0 . 0 . 0 . 0 . 0 . 0 .]; 


Mach_CD=[0 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0]; 

aoa_CD=[-180 -178 -176 -174 -172 -170 -168 -166 -164 -162 -160 -158 -135 -90 -60 -45 -30 -22 - 
20-18-16-14-12-10-8 ... 

-6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 30 45 60 90 135 158 160 162 164 166 168 170 172 


174 176 178 180]; 


cd_tab=[. 

02 .02 .02 .02 

.02 

.02 

.02 

.02 

.02; 

.03 

.03 

.03 

.03 .03 

.03 

.03 

.03 

.03; 


.05 

.05 

.05 

.05 .05 

.05 

.05 

.05 

.05; 


.08 

.08 

.08 

.08 .08 

.08 

.08 

.08 

.08; 


.11 

.11 

.11 

.11 .11 

.11 

.11 

.11 

.11; 


.14 

.14 

.14 

.14 .14 

.14 

.14 

.14 

.14; 


.185 

.185 

.185 

.185 .185 . 

185 

185 

.185 

.185; 

.235 

.235 

.235 

.235 .235 . 

235 

235 

.235 

.235; 

.25 

.25 

.25 

25 .25 

.25 

.25 

.25 

.25; 


.265 

.265 

.265 

.265 .265 . 

265 

265 

.265 

.265; 

.295 

.295 

.295 

.295 .295 . 

295 

295 

.295 

.295; 

.36 

.36 

.36 

36 .36 

.36 

.36 

.36 

.36; 



1.1945 1.1945 1.1945 1.1945 1.1945 1.1945 1.1945 1.1945 1.1945; 
2.022 2.022 2.022 2.022 2.022 2.022 2.022 2.022 2.022; 

1.662 1.662 1.662 1.662 1.662 1.662 1.662 1.662 1.662; 

1.194 1.194 1.194 1.194 1.194 1.194 1.194 1.194 1.194; 

.6 .6 .6 .6 .6 .6 .6 .6 .6; 

.3438 .3438 .3885 .4065 .414 .458 .479 .497 .514; 

.2723 .2723 .3281 .3506 .36 .415 .441 .463 .486; 

.2007 .2007 .2678 .2948 .3267 .372 .403 .43 .457; 

.1292 .1292 .2073 .2388 .2887 .329 .3655 .397 .428; 

.0576 .0576 .147 .183 .246 .286 .3278 .363 .399; 

.0174 .0174 .0225 .12 .191 .243 .29 .33 .37; 

.008 .008 .0132 .068 .127 .177 .225 .262 .297; 

.0082 .0082 .0095 .0206 .07 .113 .16 .203 .248; 

.0079 .0079 .0085 .0097 .026 .06 .1 .149 .202; 

.0075 .0075 .008 .008 .0125 .03 .065 .115 .152; 

.0075 .0075 .008 .0075 .0085 .012 .028 .066 .117; 

.0075 .0075 .008 .0075 .008 .008 .017 .05 .09; 

.008 .008 .0082 .0075 .0075 .0105 .04 .08 .1175; 

.0085 .0085 .0085 .008 .011 .036 .09 .12 .1525; 

.009 .009 .0105 .011 .029 .081 .128 .167 .203; 

.011 .011 .014 .026 .0743 .126 .17 .21 .249; 

.017 .017 .021 .08 .1247 .162 .225 .262 .298; 

.026 .026 .0935 .153 .18 .238 .285 .3225 .363; 

.145 .145 .1635 .2121 .246 .284 .326 .357 .393; 

.2147 .2147 .2259 .2643 .2887 .329 .3655 .391 .423; 

.274 .274 .2836 .3166 .3267 .327 .403 .43 .457; 

.3333 .3333 .3414 .3688 .36 .415 .441 .463 .486; 

.2927 .2927 .3991 .421 .414 .458 .479 .497 .514; 
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.6 .6 .6 .6 .6 .6 .6 .6 .6; 

1.194 1.194 1.194 1.194 1.194 1.194 1.194 1.194 1.194; 

1.662 1.662 1.662 1.662 1.662 1.662 1.662 1.662 1.662; 

2.022 2.022 2.022 2.022 2.022 2.022 2.022 2.022 2.022; 

1.1945 1.1945 1.1945 1.1945 1.1945 1.1945 1.1945 1.1945 1.1945 


36 

.36 

.36 

36 

36 .36 

.36 

.36 

.36; 


295 

.295 

.295 

.295 

.295 

295 

.295 

.295 

.295; 

265 

.265 

.265 

.265 

.265 

265 

.265 

.265 

.265; 

25 

.25 

.25 

25 

25 .25 

.25 

.25 

.25; 


235 

.235 

.235 

.235 

.235 

235 

.235 

.235 

.235; 

185 

.185 

.185 

.185 

.185 

185 

.185 

.185 

.185; 

14 

.14 

.14 

14 

14 .14 

.14 

.14 

.14; 


11 

.11 

.11 

11 

11 .11 

.11 

.11 

.11; 


08 

.08 

.08 

08 

08 .08 

.08 

.08 

.08; 


05 

.05 

.05 

05 

05 .05 

.05 

.05 

.05; 


03 

.03 

.03 

03 

03 .03 

.03 

.03 

.03; 


02 

.02 

.02 

02 

02 .02 

.02 

.02 

.02]; 



a=alpha.*180/pi; 
for i = l;leneth(a) 
ifa(j)<-180 
a(j) = a(i) + 360; 
elseif a(j) > 180 
aG) = aO) - 360; 
end 
end 

mach = abs(mach); 

CL = interp2(Mach_CL,aoa_CL,cl_tab,mach,a); 
CD = interp2(Mach_CD,aoa_CD,cd_tab,mach,a); 
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APPENDIX M. VR12CLCD.M 


function [CL,CD]=vrl2clcd(alpha) 

% vrl2clcd calculates CL and CD for the VR-12 airfoil 
% given angle of attack (alpha) in radians 
% [CL,CD]=vrl2clcd(alpha) 

CL=zeros(size(alpha)); 

CD=zeros(size(alpha)); 
a=alpha* 180/pi; 

chk=(a>=20 & a<=180); 

CL=CL+chk.*(1.1733-0.018879*a+1.5752e-3*a.^2-3.1925e-5*a.''3+2.0949e-7*a.M-4.3807e- 

10*a.^5); 


chk=(a>=-180 & a<=-50); 

CL=CL+chk.*(-4.6183-0.1923*a-3.5554e-3*a.^2-3.3273e-5*a.^3-1.4528e-7*a.M-2.3003e- 

10*a.^5); 

chk=(a>-50 & a <-30); 

CL=CL+chk.*(-0.22114+0.020857*a+2.8571e-4*a.^2); 


chk=(a>=-30 & a <=-10); 

CL=CL+chk.*(-l.ll-0.12383*a-0.01515*a.^2-6.8667e-4*a.^3-le-5*a.M); 
chk=(a<20 & a>-10); 

CL=CL+chk.*(0.11976+0.12341*a+5.5841e-4*a.''2-2.0652e-4*a.^3); 
chk=(a>=17 &a<=180); 

CD=CD+chk.*(-0.26376+0.017917*a+6.9927e-4*a.''2-9.1137e-6*a.''3+2.6277e-8*a.M); 


chk=(a>=-180&a<=-10); 

CD=CD+chk.*(-0.17486-0.034463*a-1.0233e-4*a.^2-2.8958e-6*a.^3-4.6577e-8*a.M-1.5557e- 

10*a.''5); 

chk=(a>-10 & a<=0); 

CD=CD+chk.*(9.8678e-3+3.4934e-3*a+1.4844e-3*a.''2-1.3564e-4*a.^3-1.0936e-5*a.M); 
chk=(a>0 & a<=15); 

CD=CD+chk.*(9.8e-3+7.0457e-4*a+5.6104e-5*a.''2-4.1151e-5*a.''3+3.8695e-6*a.M); 
chk=(a>15 & a<17); 

CD=CD+chk.*(-1.33+1.325e-1 *a-2.5e-3*a.^2); 
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